An Aridity Index‐Based Formulation of Streamflow Components

Direct runoff and baseflow are the two primary components of total streamflow, and their accurate estimation is indispensable for a variety of hydrologic applications. While direct runoff is the quick response stemming from surface and shallow subsurface flow paths and is often associated with floods, baseflow represents the groundwater contribution from stored sources (e.g., groundwater) to streams and is crucial for environmental flow regulations, and water supply, among others. L'vovich (1979) proposed a two‐step water balance partitioning, where precipitation is divided into direct runoff and catchment wetting, followed by the disaggregation of the latter into baseflow and evapotranspiration. Here, we investigate the role of the aridity index (ratio between mean‐annual potential evapotranspiration and precipitation) in controlling the long‐term (mean‐annual) fluxes of direct runoff and baseflow. We present an analytical solution beginning with similar assumptions as proposed by Budyko (1974), leading to two complementary expressions for the two fluxes. The aridity index explained 77% and 89% of variability in direct runoff and baseflow from 378 catchments within the continental United States, while our formulations were able to reproduce the patterns of water balance partitioning proposed by L'vovich (1979) at the mean‐annual timescale. Our approach can be used to further understand how climate and landscape controls the terrestrial water balance at mean‐annual timescales, while also representing a step toward the prediction of baseflow and direct runoff at ungauged basins.


Background
Direct runoff (Q D ) and baseflow (Q B ) are the two main components of total streamflow (Q). Q D represents the fast response of the catchment to rainfall, being mainly generated by infiltration (Horton, 1933) or saturation-excess (Dunne & Black, 1970) runoff, and subsurface stormflow (Dunne, 1978;Hewlett & Hibbert, 1967). On the other hand, Q B denotes the slow response of the catchment, which is the portion of water that did not leave the catchment as direct runoff or evapotranspiration, but is stored within the system (L'vovich, 1979;Miller et al., 2016;Ponce & Shetty, 1995a). Understanding the controls on these individual responses is of great importance to advance hydrologic sciences. For example, Q D is associated with floods (Blöschl et al., 2019) and soil erosion (Morgan & Nearing, 2011), while Q B is crucial for environmental flow regulation, and water supply (Graaf et al., 2019;Miller et al., 2016).
The need to provide freshwater for the growing populations combined with the effects of climate change has motivated studies on the prediction of global-scale groundwater recharge (Döll & Fiedler, 2008) and the assessment of climate change impacts on groundwater systems (Green et al., 2011;Meixner et al., 2016;Taylor et al., 2013), more specifically on groundwater recharge (Mohan et al., 2018;Smerdon, 2017). Catchment scale formulations of baseflow can provide additional insights on those issues as baseflow can represent groundwater recharge at sufficiently long timescales: The portion of water that is not lost as direct runoff or evaporation will eventually recharge a (shallow) aquifer and exit the catchment as baseflow. This simple yet practical approach was found to be useful in several studies on groundwater recharge estimation (Meyboom, 1961;Nathan & McMahon, 1990;Walker et al., 2019;Wittenberg & Sivapalan, 1999).
While the estimation of direct runoff has a long tradition in hydrologic sciences, most notably with the curve number (SCS-CN) method (NRCS, 2004) for event-based direct runoff estimation, very few studies investigated the controls on the mean-annual direct runoff (Q D ), whereas the controls and mechanisms of baseflow generation are still not fully understood. One of the earlier studies in which both Q D and Q B have been explicitly considered into a water balance framework was presented in the interannual water balance formulation presented in L'vovich (1979). In that study, the author observed a similar behavior between the annual precipitation and direct runoff, and between-catchment wetting (amount of precipitation not leaving through direct runoff) and baseflow. L'vovich (1979) noted that such patterns were consistent with geographic location but did not attribute physically meaningful parameters to account for the observed differences. This framework was not investigated any further until Shetty (1995a, 1995b) developed a mathematical framework to account for the observed patterns observed at the interannual timescale. Although they could eventually reproduce the patterns found in L'vovich's (1979), a physical linkage between climatic and landscape properties and the assigned parameters was not attempted. More recently, Sivapalan et al. (2011) investigated the L'vovich framework through the mathematical formulation of Ponce and Shetty (1995a) for a group of 431 catchments within the conterminous US and were able to shed light on the physical meaning of the assigned parameters by analyzing their spatial distributions across different climates. Their study, however, does not establish physical linkages between the Ponce and Shetty (1995a) model parameters and climatic and/or landscape properties. A similar approach as in Sivapalan et al. (2011) was followed by Gnann et al. (2019), who also used the Shetty (1995a, 1995b) model to explain the baseflow dissimilarities between catchments within the United States and the United Kingdom.
The climatic controls on the long-term (mean-annual) streamflow (QÞ have received much attention in the hydrologic literature and is commonly explained through the Budyko framework (Budyko, 1974). In this framework, the climate is represented by the aridity index (ϕ), which is the ratio between mean-annual potential evapotranspiration (PET) and mean-annual precipitation P À Á . The Budyko framework has been widely used for global assessments of the impacts of climate change on streamflow through differentiation of the Budyko equation (or some parametric version of it) with respect to its controlling variables (Arora, 2002;Berghuijs et al., 2017;Dooge et al., 1999;Renner et al., 2012;Roderick et al., 2014) while also being used for prediction at ungauged sites (Blöschl et al., 2013). Additionally, several studies have used this framework to draw inferences on catchment behavior at mean-annual timescales by analyzing how factors other than ϕ can explain observed departures from observed data against the Budyko equation (Berghuijs et al., 2014;Donohue et al., 2007;Milly, 1994). Wang and Wu (2013) have shown that the baseflow fraction (the ratio between baseflow and precipitation at the mean-annual scale, i.e., Q B =P) and runoff coefficient (Q=P) follow similar behaviors when plotted against ϕ for 185 perennial catchments located within the continental United States Wang and Wu (2013) derived an equation to express this relationship quantitatively, providing an interesting way to study the controls of climate on Q B . However, the equation presented in Wang and Wu (2013) have not been thoroughly analyzed for the limiting case of ϕ → 0, as it will be seen in section 1.4. More recently, Gnann et al. (2019) conducted a study to understand whether ϕ can be used to predict the baseflow fraction. Their study was based on the analysis of several hundreds of catchments from the continental United States and the United Kingdom, and their results suggest that ϕ alone cannot be used for such task. Following that, the authors parameterized the Ponce and Shetty (1995a) model to investigate the controls on baseflow generation within their group of catchments.
In this study, we discuss the role of ϕ on the long-term (mean-annual) direct runoff and baseflow (Q D and Q B Þ over a wide range of ϕ and then propose an analytical derivation that leads to two expressions for Q D and Q B as functions of ϕ. We further investigate the long-term ϕ-based formulations as potential solutions for the water balance framework proposed by L'vovich (1979) when expressed in its long-term form. The paper is organized as follows. Sections 1.2 and 1.3 provide a revision of the long-term water balance proposed by Budyko (1974) and the interannual water balance proposed by L'vovich (1979), where we discuss in detail both frameworks and establish a common nomenclature. In section 1.4, we present an approach for the analytical derivation of predictive equations for Q D and Q B based on the decomposition of Q into its complementary fluxes under similar assumptions as in Budyko (1974). This approach also provides a solution for the L'vovich (1979) framework at the mean-annual timescale. Section 2 presents the data set and methods used in this study. Following that, in section 3, we evaluate the proposed derivation to fit predictive equations for 10.1029/2020WR027123 378 catchments within the conterminous United States and test the predictive capacity of the derived equations for Q D , Q B , and Q while also assessing their ability to reproduce the spatial (between catchments) patterns arising from the two-step water balance proposed by L'vovich (1979).

The Budyko Framework
Under the assumption of negligible changes in storage over sufficiently long timescales, the water balance can be written as the partitioning of P into Q and E: Several studies have been carried out in the past with a goal of understanding the overall controls on this simple partitioning. The water balance framework proposed by Budyko (1974) and others (Mezentsev, 1955;Ol'dekop, 1911;Pike, 1964;Schreiber, 1904;Turc, 1954) is still widely used, and its success lies in the observation that the fraction of precipitation that becomes evaporation (evaporation fraction, E=P) is largely controlled by ϕ If we normalize the terms of Equation 1 by P and assume this control to be translated into a functional relationship, we can write where f E is the function that relates E=P to ϕ. Equation 2 is also known as the Budyko hypothesis. When such relationship is defined, a simple formulation for the prediction of long-term streamflow can be derived as follows: In developing a functional form for f E (ϕ), Budyko (1974) observed that the following conditions must be met: ϕ→∞; E=P→1; and Q=P→0 (4) ϕ→0; E=P→0; and Q=P→1 meaning that when aridity index reaches infinity (i.e., by having PET many times greater than P), the fraction of precipitation that becomes streamflow reaches 0, while for the situation where aridity reaches 0 (i.e., by having P many times greater than PET ), this value should reach one. Several functional forms satisfying the above constraints have been proposed in the literature (Mezentsev, 1955;Ol'dekop, 1911;Turc, 1954), the most common being the one proposed by Budyko (1974): Figure 1a shows this formulation for Q=P:

The L'vovich Framework
In the framework proposed by L'vovich (1979), the annual water balance partitioning is taken as a two-step process. In the first step, the annual precipitation (p) is partitioned into annual direct runoff (q D ) and annual catchment wetting (ω): whereas ω is further partitioned into annual baseflow (q B ) and annual evapotranspiration e: Equations 6 and 7 can be combined as 10.1029/2020WR027123

Water Resources Research
Upon observing the distinct patterns between the individual components of the two-step water partitioning in catchments across different geographical locations, L'vovich (1979) developed relationships, as shown in Figures 1b-1e, where the blue lines represent the general shape of the curves in that study. It can be seen that ω and q D respond differently to the increase in p: On one hand, ω increases almost linearly at low values of p, eventually reaching a maximum ( Figure 1b). On the other hand, q D is initially 0 at the low values of p and will only be observed when a threshold value of p is exceeded. After that, q D increases rapidly with p ( Figure 1d). A similar behavior can be found for the partitioning of ω between q B and e in Figures 1c and 1e.

Climate-Based Formulations for Budyko and L'vovich Frameworks at the Mean-Annual Timescale
From Equation 2, the runoff coefficient (Q=P) can be derived as By using Equation 8 written in terms of mean-annual fluxes, the water balance can be written as Combining Equations 9 and 10, we get which demonstrates the connection between the aridity index and the two complementary partitioning indices arising from the L'vovich (1979) formulation at the mean-annual timescale. If both left-hand terms can be written as a function of ϕ, we have Following the same reasoning as Budyko (1974), we can apply limiting conditions to the relationships in Equation 14 and 15. With values of ϕ approaching infinity, Q=P will tend to 0, which leads to both Q D =P and Q B =P approaching 0 as well with ϕ→∞; As ϕ reaches 0, Q=P will be one. In this way, some combination of Q D =P and Q B =P must occur such that their sum is equal to one. As an initial assumption (which has to be further assessed with observations), one can argue that maximum values of Q D =P and Q B =P occur at such limiting conditions: Figure 2. Location of streamflow gages for the 378 catchments with complete records used in our study, color coded by the assigned value of ϕ. and; If such relationships are established, one can rewrite the variables from the L'vovich formulation as where W represents the mean-annual catchment wetting.

Catchments and Hydrological Data
The catchments selected for this study are part of the CAMELS data set (Addor et al., 2017), which is available online and contains daily time series of streamflow and precipitation, among other variables. The meteorological variables used within the CAMELS data set are described in Newman et al. (2015). The period of analysis was 30 hydrologic years (1 October through 30 September, between 1983 and 2013). The following criteria were used for removing catchments of our analysis: (i) catchments with more than 1% missing values of streamflow, (ii) catchments receiving negative values annual evaporation (E < 0), (iii) catchments with fraction of precipitation falling as snow higher than 20% (Gnann et al., 2019), and (iv) catchments with area smaller than 20 km 2 . The resulting subset comprised data for 378 catchments (Figure 2).
Daily streamflow time series were separated into daily values of direct runoff and baseflow with a oneparameter, recursive low-pass filter (Lyne & Hollick, 1979). An inherent uncertainty of our study arises from this step, as the streamflow components obtained based on hydrograph separation cannot be attributed to specific processes. Therefore, the terms direct runoff and baseflow used here are more representative of both fast and slow dynamics capture by the digital filtering method. The filtering approach used here consisted of running the filter once, forward in time (Sivapalan et al. (2011); Wang & Wu., 2013), while the filter parameter was set to 0.925 for all catchments to assure consistency. We also adopted another approach, consisting of running the filter 3 times (forward-backwards-forward) (Gnann et al., 2019), which yielded similar results, which can be seen in the supporting information.
Finally, we calculated Q, Q D , W , Q B , and E using Equation 1 and 6-8. Additionally, daily precipitation values where aggregated into mean-annual P for each of the selected catchments.

Aridity Index Calculation
We used the Reference-crop Penman-Monteith formulation for calculating daily values of PET (in mm) as where Rn is the net radiation at the surface (MJ · m −2 · day −1 ) is the heat flux into the subsurface in (MJ · m −2 · day −1 ), e and e S are respectively the actual and saturated vapor pressure (kPa · K −1 ), u is the wind speed at 2 m (m · s −1 ), T is the air temperature at 2 m (K), Δ is the slope of the relationship between saturation vapor pressure and temperature (kPa · K −1 ), and γ is the psychrometric constant (kPa · K −1 ). Rn was calculated as where Rs is the incoming solar radiation (MJ · m −2 · day −1 ), α is the albedo of the reference surface (α ¼ 0.23), and Rnl is the net longwave radiation (MJ · m −2 · day −1 ). Briefly, we computed 10.1029/2020WR027123

Water Resources Research
Equations 21 and 22 based on the procedure described in Zotarelli et al. (2009). Daily Rs and u values were extracted for the CAMELS catchments from the gridMET data set (Abatzoglou, 2013), whereas the other necessary variables used here directly provided within the CAMELS data set. Following that, we aggregated the PET into mean-annual values and computed ϕ.
In our study, we chose to estimate PET rather than use the estimates provided with the CAMELS data set, which are first presented in Newman et al. (2015). This was done since Newman et al. (2015) used a form of the Priestley and Taylor (1971) equation in which one of its variables was used as a calibrated parameter within a hydrologic model. The estimation of this parameter, therefore, will be largely affected by errors and uncertainties in model inputs, parameters, structures, among others. Furthermore, the PET estimates of Newman et al. (2015) are not reproducible, making the verification of our results impossible for other regions of the globe. We believe that a parsimonious representation of PET, and consequently ϕ, is crucial for our analysis, since the estimated value of ϕ determines the location of a catchment along the x axis of the Budyko plot, which can potentially impact our results. Figure 3 shows the mean-annual components of the L'vovich (1979) water balance. We note a clear pattern of increment in both W and Q D with P (subplots a and d), while the concavity of the relationships is similar to what can be seen in Figures 3b-3e. A threshold value for the initiation of Q D , as suggested by L'vovich (1979) at interannual scales, is also seen at our mean-annual plots (Figure 1b), while the existence of a limiting upper value of W was not found at the mean-annual scale for the selected CAMELS catchments.

Analysis of Observed Mean-Annual L'vovich Water Balance Variables
The between-catchment patterns of the second partitioning are shown in Figures 3b and 3d. We find an increasing trend in E with W , but with higher variability. This trend, however, is not preserved for values Figure 3. Location of the selected CAMELS catchments (n ¼ 378) with respect to the L'vovich water balance, at the mean-annual time scale. While P appears to exert strong control on both W and Q D (a and c, respectively), much higher variability is present in the patterns between W and E (b) and W and Q B (d). Interestingly, the suggested limits for W and E (Figure 1, B.1 and B.2) are not seen at the mean-annual scale. Note that E is maximum for a W value close to 1,400 mm.

Water Resources Research
of W around 1,400 mm and higher. In fact, not only E is highly scattered around W ∼ 1; 400 mm, it also seems to decrease, at least for its upper limit. The patterns of Q B with W are, however, more consistent, with increasing values of W leading to an increase in Q B in a similar fashion as in Figure 3e. A threshold value for the initiation of Q B can also be seen.

Functional Forms of f D and f B
These observed patterns of Q D =P and Q B =P are shown in Figure 4, where a similarity in the patterns of both ratios as a function of ϕ is apparent. Additionally, a wide range of ϕ values is observed, ranging from 0.3 to 7.7. For simplicity, we used the functional form of an exponential decay to estimate f D and f B as where a, b, c, and d are shape parameters, while δ D and δ B are shifting coefficients necessary to satisfy conditions of Equation 17: Leading to The following split sample procedure was applied to obtain best parameter values and their uncertainties: (i) A calibration subset containing half of the total sample size (n ¼ 184) is randomly picked and fitted through a Levenberg-Marquardt nonlinear least squares algorithm (Seber & Wild, 2003), yielding estimates of a, b, c, and d. The procedure is repeated 100 times. (iii) Mean and standard deviation of the coefficient of determination (R 2 ) between predicted and observed fluxes are calculated for the validation subset, as well as the mean and standard deviation of the fitted parameters. (iv) The procedure was repeated for varying values of Q D P ! max , while its final value was chosen to be the one who yielded the best combined performances for both Q D and Q B .

Resulting Equations for f D and f B and Their Predictive Performances
The fitting procedure described above allowed us to write the final functions as f D ϕ ð Þ ¼ exp −ϕ 0:77 þ ln 0:4 ð Þ The mean values and coefficient of variation (C.V.) of the fitted exponents for f D and f B are shown in the left column of Table 1, whereas Figure 5 shows a plot of Equations 27 and 28 against observed values (subplots a and b), when mean values of fitted exponents are used. Additionally, f R as predicted by the Budyko (1974) through Equation 5, is also shown as a blue line in subplot (c). The suggested functional form for f D and f B was able to reproduce the observed patterns quite well, although it is worth noting the higher scatter , when compared to the scatter of the resulting curve f R and observed Q P .
It is also interesting to note that f R (ϕ) follows a very similar trajectory as the Budyko curve, with some noticeable differences in horizontal position of the curve against the data. Such differences are somehow expected since our data set was not the same as the one used in the original work of Budyko (1974). More importantly, it is not the objective of this work to provide an exact derivation of Equation 5.
The predictive capabilities of the derived equations can be seen at right column of Table 1, as the mean and standard deviation of the R 2 calculated using the derived equations for the validation subset of each resampling round (see section 3.2). These results provide an insight on the role of ϕ as the main control of both Q D , W, and Q B , as ϕ is able to explain 77%, 96%, and 89% of the between-catchment variability of these storage and flux terms. Moreover, low C.V. of the fitted exponents indicate a robust fit of the selected functional form. Figure 6 shows how the fitted curves (based on mean values of the fitted parameters) performed against all the observed values. A simple linear regression model (black dashed line) was fitted between observed and the predicted values discussed above and a 1:1 line (red dashed line) was added for reference. The slope of the regression line suggested a low bias, around 7% for all fluxes.
We further explore the performance of the predicted relationships for the components of the L'vovich water balance at the mean-annual timescale (Figure 7). The black circles represent the observed values for each catchment whereas the red circles are for the predicted values. Subplots (a) and (c) show a remarkable similarity between the observed and predicted variables of the first partitioning, confirming that the aridity index is able to explain great portion of the between-catchment variability, while also reproducing the trajectories of the progression between the competitions of Q D and W. The observed patterns for the second partitioning are shown in subplots (b) and (d), from where similar conclusions can be drawn. The somehow scattered pattern in the relationship between W and E seemed to be well reproduced by the aridity index formulations (subplot b). It is worth noting that the decrease in E for W ∼ 1,400 mm and higher is also reproduced, reinforcing the hypothesis of such decrease to be a function of climate. Additionally, the patterns in Q B as a function of W are also well reproduced both in magnitude of its variability as well as in the general shape of its increasing trajectory.

Discussion
Our findings highlight the role of climate on both direct runoff and baseflow at the mean-annual timescale. The proposed approach differs from the one presented by Sivapalan et al. (2011) and Gnann et al. (2019), where a conceptual model for the partitioning was implemented. The importance of our results lies in the fact that there are very few catchment scale formulations for the prediction of direct runoff and baseflow at the mean-annual timescale. Prediction of these fluxes at the catchment scale is of paramount importance. Furthermore, our approach was able to encapsulate the impact of aridity index on both direct runoff and baseflow contributions to the total streamflow, which has not been studied previously.  While enough evidence and mathematical frameworks have been thoroughly documented to explain and quantify different mechanisms in generating direct runoff at event timescales, to our knowledge, no study has shown how the long-term climate affects this hydrologic response. Moreover, our findings indicating that 77% of the variability of Q D between the selected catchment is explained by aridity index, providing an indication on how climate at short and long timescales are correlated. Direct runoff is traditionally formulated at the event scale as a function of storm characteristics, land surface properties, and antecedent moisture conditions, such as in the Curve Number method (NRCS, 2004). Guswa et al. (2018) developed an approach based on the SCS-CN method with exponential distribution of rainfall depths to predict monthly and annual values of direct runoff from 97 catchments in the conterminous United States. However, their approach requires the estimation of the CN, which is a landscape-based parameter. While their method was able to predict with high accuracy the mean-annual direct runoff when using CN obtained on rainfall-streamflow analysis of the intended catchments, it performed poorly while using the readily available tabulated values of CN. We believe the aridity index to be an emerging descriptor of the interplay between event rainfall (represented in P) and antecedent moisture conditions, which is among other factors controlled by PET.
We have shown that the aridity index alone explains 89% of the variability of mean-annual baseflow. We confirmed the results shown in Wang and Wu (2013) regarding the similarity between f R and f B . However, their formulation provides a solution for f B only and does not consider the limiting conditions for both baseflow and direct runoff coefficients at ϕ → 0. It is easy to observe in the equation from Wang and Wu (2013) that Q B =P→Q=P at ϕ → 0; that is, their solution assumes that at ϕ → 0, all precipitation becomes baseflow and no direct runoff is produced, which does not corroborate with the findings from this study (see Figure 4a).
Our results present an alternative approach to the study of Gnann et al. (2019), whose conclusion was that "baseflow fraction cannot be modelled as a function of an aridity index alone". Their study showed that catchments located within their U.K. subgroup having very low values of ϕ demonstrated no trends of increasing Q B =P with the decrease in ϕ, while a visual inspection of Figure 1b in Gnann et al. (2019) might even suggest the contrary, that is, that baseflow fraction decreases with decreasing aridity. Such behavior was explained by the existence of an upper limit for catchment storage, called the wetting potential (W P ), following the terminology of the Ponce and Shetty (1995a) model used in their work. While physically reasonable, such storage limit did not seem to play a major role in the catchments used our study, and we could not infer a negative trend in baseflow fraction with decreasing aridity (Figures 4 and 5). On the other hand, our analysis provides a simple explanation for the baseflow fraction not increasing indefinitely with We believe the major difference between both studies to be the range of ϕ used in the analysis: Our data set does not have catchments with very low ϕ values, which in Gnann et al. (2019) were shown to have been assigned lower W P values. While our study does not use the Ponce and Shetty (1995a) model, additional evidence that such low values of W P are unlikely to be found in the U.S. catchments was provided by Sivapalan et al. (2011), for which W P of very high magnitudes were also computed for the MOPEX (Duan et al., 2006) catchments as seen in Table 1 of that study.
It is possible to note from our study that the relationship between Q B =P and ϕ is subject to considerable variability, especially in the lower range of ϕ of our study. Therefore, more studies considering very humid catchments might be still needed to investigate whether the pattern seen in Gnann et al. (2019) can be seen as the expression of an even higher variability of baseflow fraction over very low aridity values.

Controls of L'vovich Water Balance Formulation
While other studies have so far used the modeling framework of Shetty (1995a, 1995b) (Gnann et al., 2019 andSivapalan et al., 2011), the relationships from Equation 18-21 were able to reproduce the overall patterns of partitioning and are based on semiempirical relationships. Even though our analysis was performed at the mean-annual scale, it suggests the use of the aridity index in interannual formulations of the water balance according to the L'vovich framework, as proposed in Sivapalan et al. (2011).

Broader Implications of the Proposed Method
Continental-and global-scale assessments of the impacts of climate change based on catchment scale formulations of the water balance have a long tradition in hydrology. Most approaches have used the Budyko framework to assess the sensitivity of changes in ϕ or its individual components on Q (Arora, 2002;Berghuijs et al., 2017;Dooge et al., 1999;Renner et al., 2012;Roderick et al., 2014). These studies apply differentiation techniques to either the Budyko equation or some other parametric version of it to derive equations relating Figure 7. Comparison between the mean-annual components of the L'vovich water balance observed at the CAMELS catchments (black) and predicted by using Equations 28 and 29 into Equations 18-21. A good agreement can be seen in (a)-(d): Both patterns of increase in W and Q D with P are reproduced (a and c), while similar patterns of variability of E as a function of W were reproduced, including the decrease of E after a threshold of W e1; 400 mm. The increase in Q B with W is also observed.

10.1029/2020WR027123
Water Resources Research the effect of changes in aridity index or its separate components (Berghuijs et al., 2017) on total streamflow. It stands to reason that our proposed approach can undergo similar procedures, yielding assessments of the effects of climate change on Q D and Q B . This might prove valuable for studies of the effects of climate on groundwater resources, as it has been previously shown in section 1.
Another interesting venue for research is the investigation of other controlling factors on Q B =P and Q D =P at the mean-annual scale. Several studies that attempted to quantify departures from the Budyko curve by use of additional landscape and other climatic factors yielded invaluable insight on how such factors affect the long-term water balance (Berghuijs et al., 2014;Donohue et al., 2007;Milly, 1994) while also improving the predictive capacity of Q D and Q B (Xu et al., 2013;Zhang et al., 2015). A similar procedure can be extended to the formulations presented here to further investigate the controls on baseflow and direct runoff coefficient. As it has been discussed here, the results from Gnann et al. (2019) point out to the importance of relating storage capacity to baseflow production, especially for very humid catchments. It would be interesting to explore how additional variables related to the catchment storage might improve the predictive capacity of the formulations shown here. Additionally, we could start by asking how the already understood climatic and landscape features that are known to provide further insight into the controls on Q are translated into the control of its complementary components Q D and Q B , which could provide a deeper insight on catchment functioning. Finally, our proposed approach may be seen as a step toward the prediction of both fluxes in ungauged basins, in which the aforementioned additional landscape and climate parameters that improve the equations shown here will become of great importance.

Conclusion
The understanding and prediction of the water balance beyond its traditional form, which considers the total streamflow as its single response can yield invaluable insight needed for different hydrologic applications. While the two-step water balance formulation proposed by L'vovich (1979) appears as a promising venue for approaching this task, the need for the understanding of its underlying climatic and landscape controls has not yet allowed for the development of robust predictive tools.
We have provided a derivation of analytical expressions of the control of ϕ on both Q D and Q B . The formulations presented here were able to explain most of the mean-annual (between-catchment) variabilities of Q D , Q B , and Q and can be seen as a valid initial step toward the prediction of these fluxes in ungauged locations. Additionally, our solution was able to reproduce the observed patterns between the components of the L'vovich (1979) water balance at the mean-annual timescale and is based on the derivation of two complementary curves that when combined provide a solution for the Q that is similar to the widely known Budyko (1974) formulation. While further investigations for assessing the ϕ-based expressions when dealing with catchments with ϕ values lower than the observed here, the large sample size (n ¼ 378) and wide range of aridity indices (from 0.3 < ϕ < 7.7) of study are encouraging for further exploring our method for regions beyond our study area.
Finally, we believe our method provides an extension for the assessments on how factors other than ϕ control the mean-annual runoff, this time allowing for the assessment of such factors on Q D and Q B and can also be used for estimating sensitivities of both components to changes in climate.