Catchment processes can amplify the effect of increasing rainfall variability

By filtering the incoming climate signal when producing streamflow, river basins can attenuate—or amplify—projected increases in rainfall variability. A common perception is that river systems dampen rainfall variability by averaging spatial and temporal variations in their watersheds. However, by analyzing 671 watersheds throughout the United States, we find that many catchments actually amplify the coefficient of variation of rainfall, and that these catchments also likely amplify changes in rainfall variability. Based on catchment-scale water balance principles, we relate that faculty to the interplay between two fundamental hydrological processes: water uptake by vegetation and the storage and subsequent release of water as discharge. By increasing plant water uptake, warmer temperatures might exacerbate the amplifying effect of catchments. More variable precipitations associated with a warmer climate are therefore expected to lead to even more variable river flows—a significant potential challenge for river transportation, ecosystem sustainability and water supply reliability.


Introduction
The temporal variability of stream flow mediates a variety of social and ecological outcomes. For example, daily flow variability determines the suitability of aquatic habitats (Fabris et al 2018), whereas variations over longer time scales affect the resilience of water supply (Vogel and Bolognese 1995), river transportation (Marengo et al 2013), local economic development (Brown and Lall 2006) and the potential for violent conflicts (Roche et al 2020). The coefficient of variation of stream flow (CV Q , defined as the ratio between the standard deviation of flow and its mean) plays a particularly important role, demarcating whether riverine processes are variance-dominated, with long periods of little to no flow interspersed with erratic bursts of high discharge, or mean-dominated, with flow rates persistently at or near their long-term mean. This distinction has implications for the form, function and resilience of river-dependent systems (Botter et al 2013).
Although driven by the variability of incoming precipitation, stream flow variability is ultimately determined by physical processes that take place throughout (and below) the land surface. Through these processes, catchments regulate stochastic weather fluctuations to sustain stream-dependent social and ecological systems, and to potentially buffer these systems against changes in these fluctuations (Chezik et al 2017, Teutschbein et al 2018. This buffering of water variability is commonly deemed an ecosystem service provided by the catchment (Guswa et al 2017), and emerges from a long-term co-evolution between the landscapes and the systems that depend on them, against the backdrop of a continually changing climate (Dietrich and Perron 2006, Sivapalan 2006, Porder 2014, Troch et al 2015, Fan et al 2019. Yet today's climate is changing at an unprecedented rate. The temporal variability of rainfall is projected to increase in most regions of the world, where increased temperatures will be associated with more intense and less frequent precipitation events (Donat et al 2016). The standard deviation of daily rainfall is projected to be more sensitive than rainfall averages to changing temperatures (4%-5% per • C versus ⩽1%-2% per • C) (Donat et al 2016, Pendergrass et al 2017. The coefficient of variation of both rainfall and stream flow has increased significantly in most regions of the United States over the last four decades (red in figure 1(a)). The ability of social and ecological systems to adapt to these changes is determined by the capacity of river catchments to buffer the effect of rapidly changing rainfall variability as they generate flow.
Using daily stream flow and precipitation data from 671 catchments located across the contiguous United States for the 1980-2015 period (see section 4), this paper highlights three important characteristics of catchments' ability to modulate climate variability. First, we show that many catchments in fact amplify the variability (as quantified by CV) of incoming precipitation (figures 1(a)-(c)). Second, these catchments might also amplify projected changes in rainfall variability, which has potentially troubling implications for stream-dependent social and ecological systems. Third, we show that both characteristics are time-scale dependent, meaning that a catchment can simultaneously dampen the CV of daily rainfall while amplifying the CV of monthly rainfall. This implies that a catchment's propensity to amplify changes in rainfall variability might depend on the time scale that is relevant for the considered application. We build on a widely used model of catchment-scale water balance dynamics (Rodriguez-Iturbe et al 1999, Porporato et al 2004, Botter et al 2007, to relate the amplifying effect of catchments to the interplay between two fundamental hydrological processes: the partitioning of rainfall into runoff and evapo-transpiration, and the retention of non-evapotranspired water as storage prior to its release as stream flow. We show that an increase in the proportion of precipitation lost to evapo-transpiration (e.g. due to higher temperatures (Berghuijs et al 2014a)) will exacerbate the amplifying effect of catchments and ultimately increase the variability of stream flow, even if rainfall variability is held constant. This suggests that projected increases in rainfall variability and evapotranspirative losses might both contribute to increasingly variable stream flows. We show that these two mechanisms had a driving influence on observed historical changes in CV Q .

Results
The ratio r = CVQ CVP between the coefficients of variation of stream flow and precipitation spans several orders of magnitude across the conterminous United States. For about a third of the 671 × 4 considered combinations of catchments and seasons, daily stream flow has a higher coefficient of variation than the incoming daily precipitation ( figure 1(b)).
However, the analysis also reveals a persistent pattern of variations in the r ratio across seasons, observation time scales and geographic regions, as shown on figure 1(c). Catchments in seasons (summer) or locations (Great Plains and Southwest) where temperature and precipitation peaks are in phase see a larger share of their precipitation 'lost' as evapotranspiration (Berghuijs et al 2014b) and tend to amplify the variability of rainfall (r > 1). In contrast, catchments where flow generation is governed by long term (seasonal) water storage, either as snow pack (Rocky Mountains and High Plains) or in the subsurface (West Coast summers), tend to dampen the variability of daily rainfall. However, a majority (83%) of catchments and seasons amplify the variability of monthly rainfall (obtained from daily observations using a 30 d moving average, figure 1(d)). Based on these results, we hypothesize that the partitioning and retention of rainfall by catchments, along with the considered observation time scale, play an important role in determining whether catchments amplify or attenuate incoming rainfall variability.
We formalize this hypothesis by deriving the r ratio of catchments based on four common, albeit extremely simplifying, assumptions about the underlying hydroclimatic processes. First, rainfall is assumed to follow a stationary marked Poisson process with exponentially distributed event depths (Rodriguez-Iturbe et al 1999). Second, instantaneous evapo-transpirative losses from the unsaturated zone are assumed proportional to its volumetric water content. Under these conditions, pulses of deep infiltration from the unsaturated zone to the saturated zone (here designated as 'recharge events') themselves follow a stationary marked Poisson process with a lower event frequency than rainfall, but with identically distributed event depths (Porporato et al 2004). Third, stream flow generation by the saturated subsurface is assumed proportional to the volume of stored water (with proportionality constant, or inverse mean hydraulic response time, k [T −1 ]), implying an exponential decrease in stream flow through time between recharge events (Botter et al 2007). Lastly, contributions to stream flow from overland runoff and lateral flows in the unsaturated zone are negligible, or effectively captured via the exponential recession model. Under these assumptions, the (squared) r ratio of catchments can be expressed as (see section 4): where f(ψ) = ψ−1+e −ψ ψ ∈ [0, 1] is a strictly increasing function. Parameter ϕ ∈ [0, 1] represents the water yield, that is, the proportion of the incoming rainfall that leaves the catchment as stream flow. Parameter ψ ⩾ 0 is the ratio between the observation time scale (e.g. monthly versus daily flows) and the mean hydraulic response time of the catchment (see section 4). Equation (1) allows isovalues of r to be mapped on the ϕ × ψ plane (figure 2(a)). Of particular interest is the isovalue line r = 1, which separates catchments that amplify (above the line on figure 2(a)) or attenuate precipitation variability. The four assumptions that underpin equation (1) are restrictive and likely fail to capture some of the processes that dominate flow generation in individual catchments. For example, the model is parametrized independently for each season, which allows it to capture seasonal changes in rainfall and temperature: a different water yield value is estimated for each season. However, carryover water storage between wet and dry seasons (both as soil moisture and groundwater) might violate the model's stationarity assumption (Müller et al 2014). In addition, the stream flow recession in many catchments might be betterapproximated as a nonlinear power-law relationship (Patnaik et al 2018), rather than the linear storagedischarge relationship assumed in equation (1). However, these processes are unlikely dominant controls on CV Q at the daily to monthly time scales that we consider, as suggested by the numerical simulations presented in the Supplementary Discussion (available online at stacks.iop.org/ERL/16/084032/mmedia).
Applied to the 671 × 4 combinations of catchments and seasons of the dataset at daily, weekly, and monthly observation time scales, the model predicted observed values of r with a mean absolute percentage error (MAPE) of 49%. As expected, predictions are substantially better in small catchments that are well aligned with the lumped nature of the model, and in regions and seasons where dominant hydrologic processes are expected to align with the theoretical assumptions of the model (figure 2(c)). For example, snow-dominated hydrology throughout the winter season in the Rocky Mountain region violates the assumption that catchment storage occurs in the subsurface and is proportional to discharge at the outlet. In many of the seasonally dry western regions, wet season onset during the Fall months of September, October, and November violate the assumption that rainfall statistics within a given season are stationary (Müller et al 2014). In the summer months of the desert Southwest, monsoonal rainfall likely arrives in short, intense bursts that trigger overland flow (Howes and Abrahams 2003), which is not strictly accounted for by the underlying runoff generation model (Botter et al 2007). Removing the 15% of catchments in region-seasons, where hydrology is likely dominated by snow (Fall and Winter in the Rockies, High Plains, Great Plains, Great Basin and Sierra Nevada) or overland flow (Summer and Fall Monsoon in the Southwest) reduces the MAPE to 33%. Notably, although simplifying assumptions on streamflow recession have little effect at the daily time scale, the theory does not accurately predict r at the weekly or monthly time scales for catchments with strongly non-linear storage discharge relationships (figure S6). We speculate that this is because longer observation time scales (and the associated time averaging) tend to shave off high flows and increase the preponderance of lower flows, for which the effect of non-linear storage-discharge on the recession limb is most salient (Karst et al 2019). Despite its shortcomings in predicting the actual value of r in some catchments, the model predicted whether catchments amplify (r > 1) or attenuate (r < 1) rainfall with an accuracy above 80% for the full sample of catchments and seasons (figure 2(b)). Model performance did not vary substantially across observation time scales (figure S6).
These results suggest that the amplifying effect of catchments is ultimately determined by a competition between two fundamental processes captured by the model, precipitation partitioning and storage retention, which appear to transcend the complex and highly local nature of stream flow generating processes in individual catchments. On the one hand, the water yield ϕ quantifies how infiltrated precipitation is partitioned between river discharge and evapotranspiration. If all other conditions remain the same, an increase in water yield increases the mean value of stream flow and therefore decreases its coefficient of variation and r ratio. On the other hand, the parameter ψ quantifies how non-evapotranspired water is retained between precipitation events and gradually released into the stream. Catchment retention affects stream flow variability differently for different observation time scales. If the observation time scale exceeds the average response time of the catchment (ψ > 1), the serial correlation introduced by catchment retention is 'integrated out' in the observed stream flow time series and will have little effect on r. In contrast, stream flow observations taken at a short enough time scale (ψ < 1) will capture the serial correlation introduced by catchment retention. These correlations attenuate the variance of stream flow and therefore decrease r.
We turn to investigating whether catchments that amplify the variability of incoming precipitation also amplify changes in that variability. This may not be the case, for instance, if changing rainfall variability affects catchments' ability to partition or retain water through changes in vegetation cover or snow pack. To estimate these effects, we use a linear regression framework, where the r ratio is taken to be a function of relative rainfall variability (CV P ) and time, and where CV P is a function of time. Taking the full timederivative of the coefficient of variation of stream flow, CV Q , then yields a linear function of ∂CVP ∂t and CV P : The first term of the derivative represents the direct (linear) response of the catchment to changing rainfall variability, as described by the ratio r. This term demonstrates that catchments that amplify rainfall variability may also amplify changes in rainfall variability. The second term represents the indirect (non linear) effect of changing rain regimes, where the r ratio is itself a function of rainfall variability. For example, more variable precipitation might affect vegetation cover in a way that decreases the proportion of rainfall lost to evapo-transpiration (Feng et al 2015). This increase in water yield would lower the r ratio (per equation (1)), in which case ∂r/∂CV P would be negative. The last terms represent the effects of an exogenous change in the r ratio that is not elicited by a change in rainfall variability. For example, warmer temperatures will cause less precipitation to fall as snow, which has been linked to lower water yields (Berghuijs et al 2014a). This would increase r (per equation (1)) irrespective of changing precipitation variability, in which case ∂r ∂t would be positive. Such changes in catchment processes will directly influence streamflow variability (third term of equation (2)) but might also feed back to affect the variability of precipitation (last term of equation (2)), for instance by affecting precipitation recycling processes (Spera et al 2016). We use linear regressions to estimate the four terms of equation (2) based on historical changes in rain and flow variability (see section 4). We focus on a subset of 563 catchments that have 35 years of continuous daily observations available, as shown in supplementary table S1. Regression results are shown in supplementary table S2 and used to construct the graphs on figure 3. Figure 3(a) shows a general historical increase in CV Q (p < 0.01) for daily, weekly and monthly observations in all seasons, except winter when CV Q decreased (p < 0.01). Changes in CV Q are dominated by exogenous changes in r ( figure 3(a), red). This suggests that relative streamflow variability might be particularly sensitive to changes in factors (such as temperature (Berghuijs et al 2014a) and land cover (Levy et al 2018)) known to strongly affect water yield and recession behavior, which our theoretical model associates with r. The indirect effect of changes in CV P (green) generally operates in the opposite direction to the direct (blue) and exogenous (red) components. This suggests that, on average across our dataset, catchments tend to 'adapt' their response to increasing rainfall variability so as to attenuate its overall effect on the variability of river flow-a negative feedback that has been associated with hydrologic resilience in past studies (Harder et al 2015). The indirect effect of CV P on changes in CV Q is also substantially smaller in magnitude than the direct effect. The ratio between the two effects is significantly smaller in absolute value than 0.55 (p < 0.01) for all seasons and observation time scales (dots on figure 3(b)). Note that catchments in different geographic regions are lumped together in the linear regressions in order for sample sizes to be sufficiently large for statistical inference. Consequently, the results in figure 3 represent average effects across catchments. To investigate regional disparities, we replicated the analysis individually for each geographic region (but lumping across seasons and observation time scales). Results in figure S4 show that the exogenous effect varies substantially across regions, as expected by the variety of factors affecting recession and water yields across regions. However, the directions and relative magnitude of the direct and indirect effects are similar to those found in the main analysis. In addition, a Monte Carlo analysis (see section 4) carried out to simulate variations of the estimated effect across catchments within the sample shows that the direct effect remains larger than the indirect effect for nearly all (>95%) simulated instances (box plots on figure 3(b)). Together, these results suggest that catchments that currently amplify rainfall variability are also likely to amplify changes in rainfall variability.

) and E[CVP]
is the observed CVP averaged across gages. Bootstrapped 95% confidence intervals are smaller than symbol sizes so not displayed. Purple boxplots represent simulated variations across catchments: α12 and α1 are drawn from independent normal distributions with mean and standard deviation given by the linear regression estimates, and CVPs are observed individually at each gage. The ratio between the average indirect and direct effects is smaller than 0.55 for all seasons and time scales (p < 0.01). The direct effect is larger than the indirect effect for >95% of the simulations. (c) Comparison between the regression coefficient α2 (black) and empirically estimated r ratios (purple), see section 4. Most regression coefficient lie within one standard deviation (error bar) of the mean empirical observation of r (empty symbol).

Conclusion
River catchments filter incoming climate signals, and therefore determine the extent to which changing rain regimes ultimately translate into changing water availability for stream-dependent social and ecological systems. Whether a catchment amplifies or attenuates changes in rainfall variability emerges from the competition between two fundamental hydrologic functions-partitioning and retention-that are active in all watersheds. The universal nature of these drivers can facilitate the assessment of climate vulnerability in datascarce basins, which is an enduring global challenge (Blöschl et al 2019). The two parameters that describe this competition-water yield and the mean catchment response time-can be directly estimated if (even sparse) rainfall and stream flow records are available (Doulatyari et al 2015, Müller andThompson 2015).
These findings have three important implications in the context of climate change. First, by decreasing the water yield (for instance by decreasing the fraction of precipitation falling as snow (Berghuijs et al 2014a)), increased temperatures will likely increase catchments' propensity to amplify the relative variability of incoming rainfall. The model associates any relative decrease in ϕ with an equivalent relative increase in r 2 , per the inverse proportional relation between ϕ and r 2 in equation (1). This implies that, by affecting water yield, increasing temperatures will increase the coefficient of variation of stream flow, even if rainfall variability remains constant. Regression results in figure 3(a) suggest that, indeed, the effects of exogenous changes in r (as might arise from changes in temperature) dominate historical changes in CV Q . Second, partitioning (ϕ) and retention (ψ) interact non-linearly to determine the r-ratio of the catchment. By capturing this relation, the model may be used to evaluate the effect of catchment alterations − land cover changes, for instance − on the amplification or attenuation of changing rainfall variability. For example, an increase in water yield from ϕ = 0.36 to ϕ = 0.4 (perhaps associated with deforestation (Levy et al 2018)) will have a disproportionately large bearing on r 2 for catchments with short response times, where a large value of ψ maps to steep isovalues of r in the ϕ × ψ plane on figure 2(a). Third, the theory elucidates the relation between observation time scales and catchments' ability to filter incoming rainfall variability. A catchment that attenuates rainfall variability for short (e.g. daily) observation time scales will nonetheless amplify rainfall variability for sufficiently long (e.g. monthly) time scales. The appropriate time scale of observation is determined by the considered application. This implies that a given catchment can both attenuate changes in rainfall variability for some applications (e.g. fish habitat driven by daily variability), while amplifying it for others (e.g. agricultural yields driven by seasonal variability). Therefore, the application context, not only the underlying hydrologic processes, determines the vulnerability of watersheds to changing rain regimes (Müller and Thompson 2019).

Data and pre-processing
Stream flow data were obtained from the Catchment Attributes and Meteorology for Large Sample Studies (CAMELS) dataset compiled by the United States Geological Survey (Addor et al 2017). We used precipitation data from the North America Land Data Assimilation data preprocessed to match the catchments of the CAMELS dataset (Addor et al 2017). The combined stream flow and precipitation dataset is publicly available at https://ral.ucar.edu/ solutions/products/camels and provides daily rainfall and stream discharge observations from 671 gaged catchments with minimal human impact in the lower 48 US states between 1980 and 2015. Observations from 16 season-catchment combinations, most of them in winter, (supplementary table S1) were removed because we were unable to identify suitable recessions to determine k (see below), likely due to snow-dominated runoff processes. Separately, we removed 108 gauges from the regression analysis due to incomplete (interrupted) time series of observations (table S1). Our sample sizes were therefore 2668 and 2252 catchment-season combinations for the validation and regression analyses, respectively. All data used in both analyses are freely available at https://curate.nd.edu/show/bc386h47534. As seen on supplementary figure S1, the dataset covers a very wide range of catchment sizes, topography, vegetation and hydroclimatic characteristics (Addor et al 2017).
Time series of daily precipitation and stream flow observations were split into four seasons according to their observation month: December to February, March to May, June to August and September to November, respectively for Winter, Spring, Summer and Fall. Moving average window of seven and thirty days were then applied to aggregate the 35 years of daily precipitation and stream flow time series into weekly (T 0 = 7) and monthly (T 0 = 30) observation time scales, respectively. The coefficients of variation (CV) of precipitation and stream flow, and their ratio r, were then computed using the aggregated time series. Note that CV Q [−] is unit-less and represents the CV of both specific (i.e. area-normalized) and total discharge. Water yields ϕ were computed for each season and catchment by taking the ratio between total seasonal stream flow volume and total seasonal precipitation volume, both taken over the whole period of data observation. Lastly, the linear recession constant k was estimated by identifying suitable recessions (at least 4 consecutive days with a decreasing, concave-up hydrograph) to be fitted with non-linear least squares, as detailed in (Dralle  et al 2017b). The recession constant was used to compute ψ = k · T 0 for each combination of catchments, seasons, and observation time scales.
We estimated temporal changes in CV P and CV Q by assuming that long-term precipitation and stream flow dynamics emerge from stationary processes that take place within multiple juxtaposed periods (Porporato et al 2006, Botter et al 2013. The 35 years of continuous daily observations available for 563 catchments of the data-set (see supplementary table S1) were split into seven non-overlapping periods of five years. To ensure that our results are not driven by the arbitrary period length, we reproduced our analyses using period lengths of two and 10 years (see supplementary figures S2 and S3). The coefficients of variation of precipitation and river flows were then computed for each period, catchment, season and (daily, weekly and monthly) observation time scales. Pairs of successive periods (t − 1 and t) were then combined to estimate temporal changes in stream flow and precipitation variability (∆CV ) across the pair.

Hydrological model
The model assumes that precipitation is a marked Poisson process, with frequency λ P [T −1 ] and exponentially distributed event volumes with (areanormalized) mean α P [L] (Rodriguez-Iturbe et al 1999). Due to the assumed independence of rain events, the CV of total rain volume, aggregated over T 0 days, is CV P = 2/(T 0 λ p ) (see Supplementary Discussion for derivation details). All incoming water is assumed to infiltrate into a subsurface unsaturated zone, meaning that canopy interception processes are embedded in α P .
Once infiltrated, water exits the unsaturated zone layer, either via evapo-transpiration (at a rate assumed proportional to current unsaturated zone water content (Porporato et al 2004)) or via percolation into the saturated zone. Such runoff generating recharge events are created when a contemporaneous rain event causes unsaturated zone moisture content to exceed an effective field-capacity, where water freely drains to the saturated zone. Under the above assumptions about the rainfall process and drainage from the unsaturated zone, it can be shown that the probability distribution of depths of recharge events is approximately equal (Botter et al 2007) to the probability distribution of rainfall depths (i.e. both distributions are exponential with mean α p ). Climate processes and water storage dynamics in the unsaturated zone therefore jointly determine the water delivered to the saturated zone, which generates discharge in the stream as the water table lowers between recharge events. These recharge events can be described as a marked Poisson process with a censored frequency λ < λ P . The censoring ratio ϕ provides a simple measure for quantifying the fraction of total rainfall that exits the catchment as stream flow over the long term (Doulatyari et al 2015), also known as the water yield: where µ p and µ Q are the long-term mean precipitated and discharged water volumes. Note that the above expression neglects contributions to total stream flow volume that are not associated with discharge generated by a subsurface saturated zone. This explains the relatively poorer performance of the model in catchments (e.g. SWS in summer, figure 2(c)) where overland flow is likely substantial. Discharge at the catchment outlet is assumed proportional to storage within the saturated zone, with a proportionality constant k [T −1 ], implying an exponential flow recession between rainfall events with characteristic drainage timescale k −1 [T]. In contrast to the rainfall and recharge events, the stream flow time series is serially correlated due to the retention and slow release of storage at a rate dependent on the hydraulic response time of the catchment k −1 (Dralle et al 2017a). This process affects the variance of cumulative stream flow as follows (see Supplementary Discussion for derivation details): Expressing mean cumulative stream flow over a period T 0 as T 0 λα P leads to an expression for the coefficient of variation of stream flow (CV Q ) and the (squared) r ratio: where ψ = kT 0 > 0 is the ratio between the observation time scale (T 0 ) and the mean response time of the catchment (k −1 ). The function f(ψ) = ψ−1+e −ψ ψ ∈ [0, 1] is strictly increasing and represents the serial correlation introduced by retention, and its attenuating effect on streamflow variability.
The model was validated by computing the absolute percentage error between the r 2 ratio predicted using equation (5) based on estimated k and ϕ, and ther 2 ratio obtained from empirical estimations of CV Q and CV P : We also evaluated the model based on the frequency of correct predictions of r < 1.

Regression analysis
We used linear regressions to analyse temporal changes in CV Q and CV P . We first assessed time trends in CV between subsequent periods. We examined the sign and statistical significance of the regression coefficient of CV Q and CV P calculated for each period, against the median year of the period. The analysis was conducted independently for each season and geographical region ( figure 1(a)), with standard errors clustered by catchment (Williams 2000).
To relate changes in CV Q to changes in CV P , equation (2) is expressed as the following linear regression: where ϵ is a random error term assumed to have a mean value of zero and be independent across catchments. Variables ∆CV Q , ∆CV P and CV P are obtained from CV Q and CV P estimated for successive 5-year periods of the observation record as described above.
Assuming that these estimates are discrete approximations of the corresponding terms in equation (2), regression coefficients can be interpreted as the exogenous (α 1 ≡ ∂r/∂t and α 0 ≡ r · ∂CV P /∂r · ∂r/∂t), direct (α 2 ≡ r) and indirect (α 12 ≡ ∂r/∂CV P ) components of ∂CVQ ∂t . Regression coefficients were estimated for each season and observation time scale using ordinary least squares, with standard errors clustered by catchment (Williams 2000).
Plugging regression estimates back into equation (6) and omitting random errors allowed us to plot the average contributions of the direct, indirect and exogenous effects on the average temporal change in CV Q across catchments, for each season and considered time scale (figure 3(a)). The relative importance of the direct and indirect effects was estimated by taking the ratio between the corresponding terms of equation (6) (figure 3(b), dots). The variability of this ratio across catchments was estimated through numerical simulations by (a) sampling CV P from the set of 563 catchment-averaged observations; (b) sampling α 2 and α 12 from independent normal distributions with mean and standard deviations given by the relevant regression estimates; and (c) computing the ratio of indirect vs. direct effects as CV P α12 α2 . Box plots on figure 3(b) represent the distribution of that ratio obtained from 1000 Monte Carlo repetitions.
We carried out two robustness checks to build confidence in our regression results. First, the analysis implies that the regression coefficient α 2 can be interpreted as the average value of r across all catchments, as seen by comparing equations (2) and (6). Figure 3(c) compares regression coefficient α 2 (black) with r ratios obtained empirically from period-of-record stream flow and precipitation observations at each catchment (purple): α 2 remains within one standard deviation (error bar) of the mean empirical value of r (empty symbols) for nearly all seasons and time scales. Second, we replicated the analysis using changes in CV Q and CV P estimated over different periods lengths. Periods longer (ten years) or shorter (two years) than the preferred duration (five years) lead to increases in both the uncertainty of regression estimates and in errors on r ≡ α 2 (supplementary figure S3). However, the mean contributions of the direct, indirect and exogenous effects, and the spatial distribution of historical trends, remain similar to those presented for periods of five years (supplementary figures S2 and S3). We interpret the errors and uncertainties for longer and shorter periods as likely caused by sample size limitations. Namely, there appears to be a tradeoff on the duration of periods between (a) having enough observations within each period to accurately estimate CV P and CV Q at daily to monthly time scales and (b) having enough periods represented to implement the linear regression analysis over a large enough sample size. The sensitivity analysis in supplementary figure S3 points to the chosen period length of five years as optimal in regards to that tradeoff.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// curate.nd.edu/show/bc386h47534.