the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Recession analysis revisited: impacts of climate on parameter estimation

### Elizabeth R. Jachens

### David E. Rupp

### Clément Roques

### John S. Selker

Recession analysis is a classical method in hydrology to
assess watersheds' hydrological properties by means of the receding limb of
a hydrograph, frequently expressed as the rate of change in discharge
($-\mathrm{d}Q/\mathrm{d}t$) against discharge (*Q*). This relationship is often assumed to take the
form of a power law $-\mathrm{d}Q/\mathrm{d}t=a{Q}^{b}$, where *a* and *b* are recession parameters. Recent
studies have highlighted major differences in the estimation of the
recession parameters depending on the method, casting doubt on our ability
to properly evaluate and compare hydrological properties across watersheds
based on recession analysis of $-\mathrm{d}Q/\mathrm{d}t$ vs. *Q*. This study shows that estimation based on
collective recessions as an average watershed response is strongly affected
by the distributions of event inter-arrival time, magnitudes, and antecedent
conditions, implying that the resulting recession parameters do not
represent watershed properties as much as they represent the climate. The
main outcome from this work highlights that proper evaluation of watershed
properties is only ensured by considering independent individual recession
events. While average properties can be assessed by considering the average
(or median) values of *a* and *b*, their variabilities provide critical insight
into the sensitivity of a watershed to the initial conditions involved prior
to each recharge event.

- Article
(1045 KB) -
Supplement
(202 KB) - BibTeX
- EndNote

Accurate representations of watershed-scale hydrological processes are
urgent in a global- and anthropogenic-change perspective. Streamflow
recession analysis has been routinely used for about half a century to
assess watershed properties (e.g.,
Brutsaert and Nieber, 1977; Kirchner, 2009; Mcmillan et al., 2014) and more
recently their vulnerability to climatic and anthropogenic factors (Berghuijs
et al., 2016; Brooks et al., 2015; Buttle, 2018; Fan et al., 2019).
Recession analysis is commonly done by plotting the time rate of change in
discharge $-\mathrm{d}Q/\mathrm{d}t$ vs. discharge *Q* with bi-logarithmic axes. Theory for an idealized
single aquifer predicts a power law relationship with parameters *a* and *b*
(Brutsaert and Nieber, 1977; Rupp and Selker, 2005).

However, it has long been recognized that the accuracy in the estimation of those parameters is highly sensitive to the methods used (Chen et al., 2018; Dralle et al., 2017; Roques et al., 2017; Rupp and Selker, 2006a; Santos et al., 2019; Stoelzle et al., 2013).

Two categories of parameter estimation methods are based on: (1) the
aggregation of all observations in the space of $-\mathrm{d}Q/\mathrm{d}t$ vs. *Q*, hereafter referred to as the
“point cloud”, to describe the average watershed behavior and (2) the
identification of individual recession events in the space of $-\mathrm{d}Q/\mathrm{d}t$ vs. *Q* to look at the
variability of a watershed's response. There is a long history of recession
analysis parameter estimation using the point cloud beginning with Brutsaert
and Nieber (1977), and it remains common (e.g.,
Buttle, 2018; Liu et al., 2016; Meriö et al., 2019; Ploum et al., 2019;
Sánchez-Murillo et al., 2015; Stewart, 2015; Vannier et al., 2014; Yeh
and Huang, 2019). In recent literature there has been a shift toward using
individual recessions to estimate recession parameters (Basso et
al., 2015; Karlsen et al., 2018; Roques et al., 2017), and Santos et al. (2019) go as far as to question the validity of point cloud estimation
methods.

When Brutseart and Nieber (1977) first proposed their recession analysis
method, aquifer recession behavior was identified by fitting a lower
envelope to the point cloud, thus assuming small values of $-\mathrm{d}Q/\mathrm{d}t$ for a given
*Q* represent aquifer discharge flow, and anything larger has contributions from
faster pathways such as overland flow. This lower-envelope (LE) method of
estimating recession analysis parameters was shown to be highly subject to
artifacts arising from measurement noise and recording precision (Rupp and Selker, 2006a; Troch et
al., 1993), and improvements to fitting a lower envelope have been proposed (Stoelzle et al.,
2013; Thomas et al., 2015). An alternative fitting method wherein *b* was
estimated as the best linear fit to the point cloud was introduced by Vogel
and Kroll (1992) as the central tendency (CT). The central tendency method was
adapted by Kirchner (2009) to address the undue weight of highly uncertain
extreme points. Kirchner (2009) also suggested fitting a polynomial function
to averages within bins of the cloud data. All of these point cloud fitting
approaches fundamentally treat each computation of $-\mathrm{d}Q/\mathrm{d}t$ and *Q* as reflecting a
single average underlying curve, with deviations from a single curve
effectively treated as noise. In other studies, data have been subset by
season or month (e.g., Szilagyi
et al., 2007; Thomas et al., 2015) to examine seasonal variations in the
recession characteristics with the subsets still treated to point cloud
analyses.

In contrast, the variability in watershed response to individual recharge
events can be depicted by fitting recession parameters to individual
recession events. Several authors have observed that individual recessions
had greater *b* values than the point cloud did (Biswal
and Marani, 2010; McMillan et al., 2011, 2014; Mutzner et al., 2013; Shaw
and Riha, 2012); a larger value of *b* indicates a time rate of decline that
decreases more quickly with decreasing streamflow. Consistent with these
studies, we have also observed individual recessions that have a larger *b* value
than the point cloud fit across watersheds in the Oregon Cascades. As an
example, we present in Fig. 1 recessions for Lookout Creek, Oregon, USA,
using daily discharge data (m^{3} s^{−1}) from 1949 to 2016 (USGS station no. 14161500; United States Geological Survey) (Johnson and Rothacher, 2019; obtained from USGS,
2019). In the 66 years of data presented, a total of 1309 recession events
are identified with an average of 19 events per year. It is clear that
values of *b* for individual recession events tend to be larger than *b* for the
point cloud, particularly at lower discharges. In this example, individual
event selection criteria include recessions lasting longer than 5 d,
starting 1 d after the peak to exclude the influence of overland flow and
ending at the next precipitation event, following other studies (Biswal
and Marani, 2010; Shaw and Riha, 2012). The *b* parameter estimated using point
cloud analysis (binning average method; BA) is smaller (*b*=1.5) compared to
the median of *b* values from the individual recessions (*b*=2.8 with 50 %
of individual recessions taking values from 2.0 to 4.7; see the color bar of Fig. 1). The frequency distribution of the *b* parameter from the individual
recessions is skewed right and roughly log-normal, which suggests that *b* from
the point cloud does not represent an average or “master” recession
behavior.

For a given discharge range in Fig. 1, there appears to be multiple
individual recessions with similar values of *b* that are horizontally offset,
implying a common *b* but a variable *a* value. The offset of individual recession
events suggests that antecedent conditions may be influencing the location
of the recession curves (e.g., Rupp et al.,
2009), consistent with various theoretical definitions of *a* that include the
aquifer saturated thickness at the onset of the recession as a parameter (Rupp and Selker, 2006b). Many authors
have associated the pattern of shifted individual recessions with
seasonality (Bart
and Hope, 2014; Dralle et al., 2015; Karlsen et al., 2018; McMillan et al.,
2011; Shaw and Riha, 2012; Tashie et al., 2019). Authors describe
a generally sinusoidal relationship with larger *a* values associated with
summer months (Dralle et
al., 2015; Shaw and Riha, 2012) and a weaker seasonal relationship for
values of *b* (Karlsen et al., 2018; Tashie et
al., 2019). Seasonality associated with meteorological conditions may well
be used as a predictor of *a* or *b*, but seasonality alone fails to address the
underlying climatic conditions that control streamflow recession. Instead of
describing the variability between events based on seasonality as a proxy,
parameter estimation should focus on antecedent and meteorological
conditions that control streamflow recession in order to form a more
comprehensive physically based understanding of recession parameters (e.g., Bart and Hope, 2014;
Karlsen et al., 2018).

This paper explores the source of the offset (ln(*a*)) and slope (*b*) on
individual recessions. Using a time series of synthetic hydrographs with
known parameters, we compare different methods for estimating the recession
analysis parameters and the sensitivity to the method on the frequency and
magnitude of events that make up the hydrograph. We are particularly
concerned with how individual recessions collectively create the emergent
point cloud and seek to describe how recession parameter estimation of the
point cloud is affected by the distribution of individual recessions.

This section presents methods for (1) the definition of three synthetic hydrographs, (2) the description of recession extraction from the hydrograph, and (3) the comparisons between four fitting methods for parameter estimation applied to a discharge time series for Lookout Creek.

## 2.1 Synthetic-hydrograph methods

This paper makes use of synthetic hydrographs to explore factors that change
*b* for individual recession events as well as the inter-arrival times of
individual events that create the point cloud. Our synthetic hydrographs are
created by defining individual recession events and stitching them together
to create a long time series. Synthetic hydrographs were chosen for this
study because each individual recession can be definitively identified, as
the characteristics are known, which is unrealistic when considering real
watersheds. Furthermore, the synthetic hydrographs can be specified to
directly compare different climatic controls without the confounding
variables traditionally associated with real watersheds. For these purposes,
the specifications of the synthetic hydrographs were chosen to explore the
effects of the magnitudes and frequency of recharge events on the recession
analysis parameters from collective vs. individual recessions.

The falling limb of the hydrograph is assumed to follow a power law following Eq. (2) (Dewandel et al., 2003; Drogue, 1972; Rupp and Woods, 2008):

where *Q* is the discharge, *Q*_{0} is the initial discharge prior recession at
*t*=0, *t* is the time in days since the recession started, *τ* is a
characteristic timescale, and *w* is the dimensionless power law decay
exponent. Equation (2) can be expressed as Eq. (1) with $a=w/\left(\mathit{\tau}{Q}_{\mathrm{0}}^{\mathrm{1}/w}\right)$ and $b=(\mathrm{1}+w)/w$.

Holding *τ* constant and varying the initial condition *Q*_{0}, results in
a hysteretic relationship of $-\mathrm{d}Q/\mathrm{d}t$ vs. *Q*, in contrast to a constant *a* value which
produces a single non-hysteric relationship. Defining *a* as a function of
initial conditions has both theoretical (e.g., Rupp and Selker, 2006b) and
empirical (e.g., Bart and Hope, 2014) support.
The constancy of *τ* is not well established, but we assume it is
constant for the scenarios examined here. Consequentially, a constant
*τ* results in a variable value for *a* that is inversely proportional to the
initial discharge. An inverse relationship is consistent with theoretical
expectations for non-linear aquifers (*b*>1), where *Q*_{0}
increases with increasing initial saturated thickness (see Figs. 2 and 3 in
Rupp et al., 2006b). Though the particular timescale is not important to our
objectives, we chose it to be 45 d. Brutsaert (2008) noted a tendency for
*τ* to be near 45 d across a large number of basins when fitting Eq. (1) with *b*=1 to point cloud data. It remains to be seen whether a
similarly narrow distribution of *τ* occurs for *b* not equal to 1.

A pulse recharge amount corresponding to a given *Q*_{0} can be calculated by
integrating Eq. (2) from *t*=0 to *t*=∞. For *w*>1 (*b*<2), the recharge volume is

where *D* is the depth of recharge and *A* is the aquifer area. For *w*≤1 (*b*≥2), integrating Eq. (2) results in an infinite
volume, so *b*>2 can only be sustained over a finite part of any
recession. Values of *b*>2 have been derived from the physical theory
for the early portion of a recession (Brutsaert and Nieber, 1977; Rupp and
Selker, 2005) or can be obtained from recession curves over a finite time
period while retaining physical realism by combining discharge from multiple
linear (*b*=1) or non-linear ($\mathrm{1}<b<\mathrm{2}$) reservoirs (e.g., McMillan et al., 2011). The effect
on *b* of combining linear reservoirs in parallel (e.g.,
Clark et al., 2009; Gao et al., 2017; Harman et al., 2009) and series (e.g.,
Rupp et al., 2009; Wang, 2011) has received much more attention.

We compared three hypothetical time series generated with different
assumptions about the distribution of the magnitudes and inter-arrival times
of recharge events and the superposition of recession events (Table 1). The
inter-arrival times are distributed log-normally (Cases 1 and 3) or
uniformly (Case 2). Event magnitudes (as defined by *Q*_{0}) are either
distributed log-normally (Cases 1 and 3) or have a constant magnitude (Case 2). Events are either independent of antecedent conditions (Case 1), or
events are superimposed on antecedent conditions (Cases 2 and 3) (Table 1
and Fig. 2).

To generate the time series for Cases 1 and 3, independent recessions were
created using a random-number generator for a log-normal distribution for
event peak magnitude and duration for a total of 10 years of time series
data. The log-normal distributions for event magnitude and duration are
chosen for the synthetic hydrographs because the distributions for Lookout
Creek are skewed right and roughly log-normal (Supplement Fig. S1), which is also consistent with
other skewed-right precipitation distributions in previous studies (Begueria et al., 2009; Selker and Haith,
1990). Recharge events were created with log-normally distributed
inter-arrival times (*μ*=2.5, *σ*=1) and event magnitudes
(*μ*=1 d, *σ*=1) where both values are normalized by
timescale and the unit hydrograph respectively, resulting in dimensionless
quantities. These values of *μ* and *σ* result in event lengths with
a mean of 20 d and an average of 18 events per year. This value was
chosen to be comparable to the 19 events per year identified in the Lookout
Creek example. The distributions of both the inter-arrival times and event
magnitudes are skewed right, representing the high frequency of smaller
events and less frequent large events. Changing *μ* and *σ* will
modify the amount of variability in individual recessions and could be
further explored with different distributions in future research regarding
the resulting variability in *b*. Case 2 assumes a constant event inter-arrival
time ($\mathit{\mu}=\mathrm{450}/\mathit{\tau}$) and magnitudes (*μ*=1). The mean
inter-arrival time of 10 d is intended to be comparable with the 19
events per year identified in the Lookout Creek example.

For Case 1, the individual recessions were combined to make a time series
such that each event was concatenated onto the last event disregarding the
antecedent flows. For Cases 2 and 3, individual recessions were superimposed
on antecedent flows, appealing to the simplest model presented by the
instantaneous unit hydrograph method (Dooge, 1973). We
acknowledge that the framework for the instantaneous unit hydrograph as
described in Dooge (1973) does not consider non-linear reservoirs, but we use
it as a simple representation to produce variability between recessions. We
discuss the implicit assumptions of this model in the Discussions and
Conclusions section. From Fig. 2, the baseflow from the first event,
*Q*_{B}, is a simple continuation of the first recession. The underlying
second event, *Q*_{C}, is defined by the second event's initial magnitude
(constant in Case 2 and randomly generated in Case 3). The resulting flow,
*Q*_{D}, is the sum of *Q*_{B} and *Q*_{C}.

As a result, Case 1 looks specifically at a time series events where the
falling limb of each event maintains the same decay constant and the effect
of having no antecedent baseflow influence on streamflow. By including
baseflow to Case 2 but maintaining equal inter-arrival times and event
magnitudes, we look specifically at the effect of antecedent conditions on
individual recessions and the point cloud. Case 3 combines the distribution
of event inter-arrival times and magnitudes of Case 1 with the baseflow
conditions of Case 2, best representing the variability and inter-arrival
times of individual recession events seen in Fig. 1 for data from Lookout
Creek. Each case will address how the controls on the hydrograph affect the
recession analysis plot and the estimates of *a* and *b*.

## 2.2 Recession extraction method

Recession extraction from observed hydrographs and the associated
sensitivities to different criteria have been explored by Dralle et al. (2017), including minimum recession length and the definition of the
beginning and the end of the event. For Lookout Creek, we used extraction
criteria similar to those of other studies (e.g.,
Chen and Krajewski, 2016; Dralle et al., 2017; Stoelzle et al., 2013) and
applied the same criteria prior to all fitting methods presented in Sect. 2.3 to isolate differences in calculated *b* values due to the fitting method only.
An individual recession event duration must be longer than 5 d. Rainfall
data can be used to identify non-interrupted recessions, but rainfall data will not available in all cases, so we rely on the hydrograph only. The
start of the recession is defined as 1 d after the discharge peak to
account for the presence of overland flow. The end of the recession occurs
at the minimum discharge prior to an increase in discharge greater than the
error associated with instrument precision for stage height of
∼0.5 cm, which translates into errors in discharge from
∼0.01–0.1 m^{3} s^{−1}, depending on the rating curve and the
discharge level (Thomas et al., 2015).

For the synthetic hydrographs used in Sect. 3.2, events of any length were included; the recession start was selected at peak discharge because overland flow was not a factor; and the end of the recession was chosen as the time immediately before the next generated discharge peak.

## 2.3 Parameter estimation methods

Four methods of estimating representative recession parameters were evaluated: lower envelope, central tendency, binning average (Kirchner, 2009), and the median individual recessions (MI) (Roques et al., 2017). Linear regression in bi-logarithmic space was used with each method for consistency across methods.

Because a change of hydraulic regime was suggestive in Fig. 1 between high-flow ranges and low-flow ranges, recession analysis parameters were
estimated for two flow ranges, early time and late time. Early time and
late time describe a theoretical transition of flow regimes between
high-flow and low-flow ranges (Brutsaert and Nieber, 1977). To reduce the
subjectivity of distinguishing between high and low flows, a breakpoint in
discharge separating high- from low-flow behavior was optimized to best
represent the analytical solution. By separating the data into two
subgroups, either smaller or larger than a defined breakpoint discharge, the
best-fit line was determined for each subgroup. The location of the
breakpoint is defined so the error between the observed ratio of *b* for the
two subgroups and the theoretical ratio (*b*=3 for early and 1.5 for late
give a ratio of 2) is minimized, theoretically defining the subgroup above
the breakpoint as early time and the subgroup below the breakpoint at
late time.

For each of the four estimated methods, parameters were estimated for the
early-time and late-time behavior separately. For the LE method, *b* was fixed
to 3 and 1.5 for early and late time, respectively, following Brutsaert and
Nieber (1977), and *a* was chosen such that 5 % of points were below the lower
envelope (Brutsaert,
2008; Troch et al., 1993; Wang, 2011). It should be noted that using these
values for *b* assumes that the groundwater discharge behaves like discharge
from a single, initially saturated, and homogenous Boussinesq aquifer. An
alternative method to fitting the lower envelope without a predefined value
of *b* was introduced by Thomas et al. (2015) using quantile regression to
estimate both *a* and *b*, but it was not used in this study. For the CT method, the
fit included all points of $-\mathrm{d}Q/\mathrm{d}t$ vs. *Q* as unweighted (Vogel and Kroll, 1992). For
the BA method, bins spanned at least 1 % of the logarithmic range, and a
line, instead of the polynomial suggested by Kirchner (2009), was fit to the
binned data. We dispensed with the inverse-variance bin weighting used by
Kirchner (2009) to account for data noise when we applied the method to the
synthetic recessions because some bins contained few points with very low
variance and therefore were weighted excessively. For the MI method,
parameters were estimated for individual recessions following Roques et al. (2017), and the medians of *a* and *b* were calculated. In all cases, the time
derivative of $-\mathrm{d}Q/\mathrm{d}t$ was computed using the exponential-time-step (ETS) method
proposed by Roques et al. (2017).

## 3.1 Parameter estimation for observed recessions (Lookout Creek)

In Fig. 3 we display the recession plot stacking all individual recessions
resulting in the formation of the point cloud. The different fitting
strategies revealed that the LE, CT, and BA methods all fit to the point
cloud and result in different values of *b* when applied to the observed daily
averaged streamflow for Lookout Creek: early-time values of *b* are over 50 %
larger for LE (fixed at 1.5) than CT and BA, and late-time values of *b* are
50 % and 25 % larger for LE than CT and BA, respectively (Fig. 3 and
Table 2). The CT and BA methods are fairly consistent with each other for
both early and late time, whereas the predefined theoretical *b* values for
the LE appear to provide poorer fits to the point cloud.

More importantly, parameter estimation differs greatly whether the point
cloud or individual recessions are used. The late-time *b* value which defines
the low-flow baseflow regime is 6 times greater for MI than CT (Table 2).
Using the MI method, the *b* value is larger than any other method for both
early and late time.

## 3.2 Synthetic-hydrograph results

Based on the similar results from BA and CT methods discussed above, and the
questionable practice of setting an early- and late-time *b* a priori as we did
in the LE method, hereafter we use the BA method to represent the point cloud
recession parameter estimation when comparing it to the MI method using
individual recessions.

The recession decay exponent *w* in Eq. (2) was set to 1.2; distinct values of
*w* were not used for early and late time. This value for *w* results in *b*=1.8
for an individual synthetic recession, which is near the reported median of
individual *b* values of 2.0 in Biswal and Marani (2010) and 2.1 in Shaw and
Riha (2012) and Roques et al. (2017), though less than the median individual
*b* of 2.8 for Lookout Creek.

The *b* values and the offset of individual recessions resulting from Eq. (1)
are both functions of *w*. A larger *b* value indicates a more stable baseflow
discharge (a slower decline rate for given discharge). For a given value of
*b* and *τ*, *a* varies inversely with ${Q}_{\mathrm{0}}^{\mathrm{1}/(b-\mathrm{1})}$.
Decreasing *w* results in larger values of *b* while also increasing the offset
between individual recessions, resulting in a larger range of *a* values and a
more scattered point cloud. In contrast, as *w* approaches infinity, the offset
is minimized as *b* goes to 1, representing an exponential falling-limb
recession in time (Rupp and Woods, 2008). In this special case,
all of the individual recessions overlap with constant *a* (i.e., there is no
offset among individual recessions lines). While *b*=1 is interpreted as a
linear reservoir according to traditional theory and is a convenience often
assumed, this result suggests that a condition where *b*=1 and *a* is a constant
is not consistent with the existence of a point cloud, except to the degree
at which observation error introduces noise into the recession hydrograph
or other pathways (e.g., overland flow) contribute to the flow in the
stream. In summary, the more linear the response is (the closer *b* is to 1),
the smaller the offset is, whereas the more non-linear the response (the larger
the *b*) is, the greater the offset and thus the more different the
parameter estimations between the point cloud and individual
recession methods will be. The three following cases using synthetic hydrographs are
intended to highlight the offset of the individual recession curves.

### 3.2.1 Case 1

Recession analysis of a hydrograph with log-normally distributed event
inter-arrival times and peak discharge with a constant falling-limb decay
constant (no baseflow represented) results in individual recession events
with the same *b*, horizontally shifted based on the initial discharge (Fig. 4). In this case, the peak flow of the event is the only source of
variability in the recession parameter *a*. The variable event magnitudes
result in individual events located over a range of ln(*Q*) values, whereas if
the same flow magnitude was preserved for each event, each individual
recession would plot on top of one another creating a single line without a
point cloud. The variable event inter-arrival times change the duration of
an event, with longer events occurring over a greater range on the *y* axis.
In this simple hydrograph, parallel individual recessions are present, and all
have *b*=1.8, as expected. The value of *b* is estimated at 1.3 considering
the point cloud, which appears to be significantly less than imposed
individual *b* value of 1.8. This underestimation results from the offset between
individual recessions based on the range of initial discharges.

To examine the sensitivity of parameter estimation to recession extraction
criteria, we evaluated how choosing the start of the recession (i.e., the
time elapsed since peak discharge) affects the value of *a* when using the
point cloud method. Whether we chose 0, 1, or 2 d following peak
discharge, *a* from the point cloud *a* was 0.17 (–), and *b* was 1.3 (–).

### 3.2.2 Case 2

The superposition of recession events accounts for the effects of antecedent
baseflow. The superposition changes the effective *w* of the falling limb of
the hydrograph as the event recession is added to the antecedent events,
resulting in a variable *b* value across the individual recessions (Fig. 5). The median
*b* value represented is 3.25 with a range of 2.56 to 3.41 (quantile range
represented in the color bar of Fig. 5). The point cloud *b* of 2.35 falls
outside of the range of *b* values for individual recessions. Superposition
results in a larger *b* value than what would arise from non-superposition. Steeper
recessions (higher *b* value) are associated with events with higher baseflow
contribution given the same addition of flow. By including antecedent-flow
conditions, neither *a* nor *b* is preserved between individual recessions.

### 3.2.3 Case 3

A hydrograph more representative of real-case conditions includes variable
inter-arrival times and event magnitudes from Case 1 and baseflow antecedent
conditions from Case 2 (Fig. 6a). These complexities result in a recession
plot where the individual recessions represent the variability in watershed
response represented by the hydrograph (Fig. 6b), where *a* and *b* are different
between individual recessions. As with Cases 1 and 2, the median individual
*b* value (3.3) is greater than the point cloud *b* value (2.0). The minimum individual *b* value is 1.9
with a maximum of 8.5, while the point cloud *b* is near the low end of the
range of individual *b* values (see the color bar of Fig. 6). The similarity of
features of Figs. 6b and 1 are noteworthy. Though many of the observed
recessions in Fig. 1 are slightly curvilinear (in the log–log space),
the synthetic recessions are power laws; in both cases there is a tendency
for recessions with lower initial discharges to have higher values of *b*, yet
still many instances of recessions have similar initial discharges but
different values of *b*.

In the 42 years since Brutsaert and Nieber (1977) proposed their recession analysis, it has provided a seemingly simple analytical method for estimating basin-scale hydrologic properties. However, recent studies have highlighted the sensitivity to estimation methods on the recession parameter values and to the resulting interpretation of average watershed behavior. This paper explores the effect of the distribution (in time and in magnitude) of individual recessions on parameter estimation and compares that to the parameter estimation from collective recessions (i.e., the point cloud). The four estimation methods considered were the lower envelope, central tendency, binning, and individual recession method. Because of the poorer apparent fit and problems pointed out from previous studies when using the lower envelope and central tendency methods, we chose to use the binning method to compare with results from the individual recessions method for a set of synthetic case studies.

We hypothesize that the climate controls the distribution of individual
recessions in bi-logarithmic plots of $-\mathrm{d}Q/\mathrm{d}t$ vs. *Q*. This distribution can be
related to the variability in recession analysis parameters. Using the three
synthetic case studies, we examine the effects of event inter-arrival time,
magnitude, and antecedent conditions on the distribution of individual
recession events that comprise the pattern of collective recessions.

We conclude that recession analysis performed on collective recessions does
not capture average watershed behavior, regardless of the fitting method.
The point cloud is an artifact of the variability of the individual
recessions, including the event inter-arrival times and distribution of
magnitudes. Individual recessions with the same *b* but different *a* values can be
produced by varying the initial discharges (Case 1); variability of *b* for
individual recessions can be produced by superimposing events on antecedent-flow conditions (Case 2); and different recession event lengths with
different *b* values can be produced by including variable event inter-arrival times
and magnitudes (Case 3).

For Case 1, the recession analysis parameter *a* is equal to $w/\left(\mathit{\tau}{Q}_{\mathrm{0}}^{\mathrm{1}/w}\right)$, and thus the intercept of the individual recession curves
will scale with *Q*_{0}. The result is a collection of individual recession
curves that are horizontally offset based on the initial discharge, producing
a smaller *b* value for the point cloud compared to the individual recessions.
Case 1 illustrates that the slope of individual recession events can be
greater than the best-fit line through the point cloud, consistent with
previous studies (Biswal
and Marani, 2010; Mutzner et al., 2013; Shaw and Riha, 2012). However, the
point cloud in Case 1 is generated by a collection of multiple individual
recessions, all with the same slope, and does not have the variability in *b* values
presented in these previous studies and shown for Lookout Creek in Fig. 1.
Cases 2 and 3 use superposition of antecedent-flow events that
consequentially changes the individual *b* values, providing a possible
explanation for the variability in *b* values for individual recessions. For
Case 2, the superposition of events takes account of antecedent conditions
which results in a distribution of individual *b* values despite the decay
exponent *w* being fixed. For Case 3, the horizontal offset of individual
recessions from Case 1 and the effects of antecedent conditions from Case 2
result in the recessions with variable individual *b* values that are
horizontally offset to create a pattern similar to that observed in a real
watershed.

While the mean *b* for individual recessions in Case 1 is a direct consequence
of the value of *w* used in Eq. (2), this is not true when the discharge from
each application of Eq. (2), which we call an event, is superimposed on
the antecedent flow, as in Cases 2 and 3. This superposition of events
results in a range of individual recession *b* as often observed in the
literature (Basso
et al., 2015; Biswal and Marani, 2010; Mcmillan et al., 2014; Mutzner et
al., 2013; Shaw and Riha, 2012); thus it appears that the straightforward
superposition of events can recreate the watershed behavior. However, there
is a key underlying assumption of this superposition that is inconsistent
with a real watershed. To help describe this inconsistency, we can compare
two distinctly different conceptual models of watersheds. The first, and very
frequently used, model is a single bucket with an outlet near the bottom.
The bucket contains a porous medium whose properties may vary with depth to
create a variety of non-linear outflows. Each new recharge event adds to the
pre-existing storage of water in the bucket. The second model is the one
used for Cases 2 and 3: each new event adds water to a new and independent
bucket, and the outflows from all buckets are aggregated. Both conceptual
models have components that are patently unrealistic when applied to natural
watersheds, but, remarkably, the latter model produces a distribution of
recession events in the space of $-\mathrm{d}Q/\mathrm{d}t$ vs. *Q* that is more like what is observed in
the Lookout Creek basin and others (Mutzner
et al., 2013; Shaw and Riha, 2012). This finding reveals key information
about the subsurface plumbing system of the basin and its dynamics that
could be explored with models with a higher degree of realism. We
acknowledge that there are other ways to create watershed memory that would
also generate variability in apparent recession parameters and would be
worthwhile to consider. For example, following previous works that have
shown that multiple linear reservoirs can generate power law recessions
(Clark et al., 2009; Harman et al., 2009), one could explore combinations of
parallel linear reservoirs with varying sizes and recession constants under
time-varying recharge. However, based on the results of Harman et al. (2009)
using periodic recharge events, it is not clear that this would lead to a
distribution of recession curves of varying *b* values like what is seen in
Fig. 1. A similar, albeit more complicated, exercise could also be done with
combinations of parallel non-linear reservoirs with distinct recession
parameters.

An additional important simplifying assumption of this study is the use of a
constant timescale *τ* for each individual event. Previous studies that
have examined timescales across basins by setting *b*=1 and estimating
*τ* from the point cloud (Brutsaert, 2008;
Lyon et al., 2015). However, given the questionable of the validity of the
point cloud estimation methods, additional studies of the variability of the
timescale among individual recession events and across the basin should be done.

We show how the point cloud pattern does not arise from watershed properties
alone. The consequence is that parameters estimated from the point cloud do
not represent watershed properties. In all three synthetic-hydrograph
representations, the median individual recession *b* is significantly greater
than *b* from the point cloud. Additionally, it is possible for the point cloud
*b* to be smaller than the minimum individual recession *b* indicating the point
cloud fit represents a behavior outside the range of watershed responses
represented by individual recession events. In contrast to the point cloud,
individual recession analysis provides insights into the average and
variability of watershed responses which is highly dependent on the memory
effect of the watershed. The variability in individual responses gives
insight into watershed complexities including heterogeneity in topography,
geology, and climate. Watersheds may present large variability in geology
and so hydrogeological conditions such as unconfined or confined aquifers,
inter-basin groundwater flows, high spatial hydraulic conductivity
variability, depth-dependent hydraulic conductivity, or large-scale
discontinuities. As a result, there are still opportunities to further
characterize the variability in watershed responses and the associated
variability in individual *b* values to improve streamflow prediction using
recession analysis.

A strength of the critical-zone community is the ability to create a global
analysis by comparing across studies (Brooks
et al., 2015; Fan et al., 2019). However, a lack of consensus for a standard
method for recession analysis procedures exists and thus inhibits recession
analysis studies from being widely compared. If streamflow analysis is to be
included in a global analysis, results need to be comparable across scales
and observatories. There is a need for a common method employed to compare
the average and variability in watershed responses. Because estimated
parameters may differ greatly by estimation method, misinterpretation of
hydrological properties and incorrect predictions within the critical zone
are possible. When using the point cloud in particular, a smaller recession
parameter *b* at a late time could be, and has been, interpreted to imply greater
basin vulnerability to drought (e.g., Berghuijs et al., 2016, 2014; Yeh and
Huang, 2019). However, a more stable baseflow is implied by the distribution of
*b* from the individual recessions and its median *b* value, which can be much larger
than what is estimated from the point cloud. We suggest that the use of
collective recession analysis should be avoided in favor of individual
recession analysis as the standard to describe the average and variability
in watershed response. The methods employed for recession analysis certainly
require more attention. Correct methods are critical to understand the
underlying hydrology and thus the interpretation of a watershed's
vulnerability to climate change.

The streamflow record for Lookout Creek is freely available from the USGS website (https://waterdata.usgs.gov/nwis/uv?site_no=14161500; USGS, 2019). The source code for the exponential-time-step method is available by request (Roques et al., 2017). Randomly generated log-normal event magnitudes and inter-arrival times presented in this paper for Cases 1 and 3 are available at: https://doi.org/10.4211/hs.e3c159631acd470cbeef5fa1abe0142e (Jachens, 2020). Respective codes can be obtained from the corresponding author.

The supplement related to this article is available online at: https://doi.org/10.5194/hess-24-1159-2020-supplement.

ERJ, DER, and JSS were involved in conceptualization. ERJ and CR developed the methodology and performed the analysis. ERJ prepared the paper with contributions from all co-authors.

The authors declare that they have no conflict of interest.

This work was supported by the National Science Foundation awards to the Center for Transformative Environmental Monitoring Programs (nos. 1551483 and 1440506).

This paper was edited by Daniel Viviroli and reviewed by Michael Stoelzle and one anonymous referee.

Bart, R. and Hope, A.: Inter-seasonal variability in baseflow recession rates: The role of aquifer antecedent storage in central California watersheds, J. Hydrol., 519, 205–213, https://doi.org/10.1016/j.jhydrol.2014.07.020, 2014.

Basso, S., Schirmer, M., and Botter, G.: On the emergence of heavy-tailed streamflow distributions, Adv. Water Resour., 82, 98–105, https://doi.org/10.1016/j.advwatres.2015.04.013, 2015.

Begueria, S., Vicente-Serrano, S. M., Lopez-Moreno, J. I., and Garcia-Ruiz, J. M.: Annual and seasonal mapping of peak intensity, magnitude and duration of extreme precipitation events across a climatic gradient, northeast Spain, Int. J. Climatol., 29, 1759–1779, https://doi.org/10.1002/joc.1808, 2009.

Berghuijs, W. R., Sivapalan, M., Woods, R. A., and Savenije, H. H. G.: Patterns of similarity of seasonal water balances: A window into streamflow variability over a range of time scales, Water Resour. Res., 50, 5638–5661, https://doi.org/10.1002/2014WR015692, 2014.

Berghuijs, W. R., Hartmann, A., and Woods, R. A.: Streamflow sensitivity to water storage changes across Europe, Geophys. Res. Lett., 43, 1980–1987, https://doi.org/10.1002/2016GL067927, 2016.

Biswal, B. and Marani, M.: Geomorphological origin of recession curves, Geophys. Res. Lett., 37, 1–5, https://doi.org/10.1029/2010GL045415, 2010.

Brooks, P., Chorover, J., Fan, Y., Godsey, S. E., Maxwell, R. M., McNamara, J., and Tague, C.: Hydrological partitioning in the critical zone: Recent advances and opportunities for developing transferable understanding of water cycle dynamics, Water Resour. Res., 51, 6973–6987, https://doi.org/10.1002/2015WR017039, 2015.

Brutsaert, W.: Long-term groundwater storage trends estimated from streamflow records: Climatic perspective, Water Resour. Res., 44, 1–7, https://doi.org/10.1029/2007WR006518, 2008.

Brutsaert, W. and Nieber, J. L.: Regionalized drought flow hydrographs from a mature glaciated plateau, Water Resour. Res., 13, 637–643, https://doi.org/10.1029/WR013i003p00637, 1977.

Buttle, J. M.: Mediating stream baseflow response to climate change: The role of basin storage, Hydrol. Process., 32, 363–378, https://doi.org/10.1002/hyp.11418, 2018.

Chen, B. and Krajewski, W.: Analysing individual recession events: sensitivity of parameter determination to the calculation procedure, Hydrolog. Sci. J., 61, 2887–2901, https://doi.org/10.1080/02626667.2016.1170940, 2016.

Chen, X., Kumar, M., Basso, S., and Marani, M.: On the effectiveness of recession analysis methods for capturing the characteristic storage-discharge relation: An intercomparison study, Hydrol. Earth Syst. Sci. Discuss., https://doi.org/10.5194/hess-2018-65, 2018.

Clark, M. P., Rupp, D. E., Woods, R. A., Tromp-van Meerveld, H. J., Peters, N. E., and Freer, J. E.: Consistency between hydrological models and field observations: Linking processes at the hillslope scale to hydrological responses at the watershed scale, Hydrol. Process., 23, 311–319, https://doi.org/10.1002/hyp.7154, 2009.

Dewandel, B., Lachassagne, P., Bakalowicz, M., Weng, P., and Al-Malki, A.: Evaluation of aquifer thickness by analysing recession hydrographs. Application to the Oman ophiolite hard-rock aquifer, J. Hydrol., 274, 248–269, https://doi.org/10.1016/S0022-1694(02)00418-3, 2003.

Dooge, J. C. I.: Linear Theory Of Hydrologic Systems – Technical Bulletin No. 1468, Agric. Res. Serv. United States Deparment Agric., Washington, D.C., 1973.

Dralle, D., Karst, N., and Thompson, S. E.: A, b careful: The challenge of scale invariance for comparative analyses in power law models of the streamflow recession, Geophys. Res. Lett., 42, 9285–9293, https://doi.org/10.1002/2015GL066007, 2015.

Dralle, D. N., Karst, N. J., Charalampous, K., Veenstra, A., and Thompson, S. E.: Event-scale power law recession analysis: quantifying methodological uncertainty, Hydrol. Earth Syst. Sci., 21, 65–81, https://doi.org/10.5194/hess-21-65-2017, 2017.

Drogue, C.: Analyse statistique des hydrogrammes de decrues des sources karstiques statistical analysis of hydrographs of karstic springs, J. Hydrol., 15, 49–68, 1972.

Fan, Y., Clark, M., Lawrence, D., Swenson, S., Band, L. E., Brantley, S., Brooks, P., Dietrich, W., Flores, A., Grant, G., Kirchner, J., Mackay, D., McDonnell, J., Milly, P., Sullivan, P., Tague, C., Ajami, H., Chaney, N., Hartmann, A., Hazenberg, P., McNamara, J., Pellet, J., Volk, J., and Yamazaki, D.: Hillslope Hydrology in Global Change Research and Earth System Modeling, Water Resour. Res., 55, 1737–1772, https://doi.org/10.1029/2018WR023903, 2019.

Gao, M., Chen, X., Liu, J., Zhang, Z., and Cheng, Q.-B.: Using Two Parallel Linear Reservoirs to Express Multiple Relations of Power-Law Recession Curves, J. Hydrol. Eng., 22, 1–12, https://doi.org/10.1061/(asce)he.1943-5584.0001518, 2017.

Harman, C. J., Sivapalan, M., and Kumar, P.: Power law catchment-scale recessions arising from heterogeneous linear small-scale dynamics, Water Resour. Res., 45, 1–13, https://doi.org/10.1029/2008WR007392, 2009.

Jachens, E. R.: Data for Cases 1 & 3 (event magnitudes and inter-arrival times), HydroShare, https://doi.org/10.4211/hs.e3c159631acd470cbeef5fa1abe0142e, 2020.

Johnson, S. and Rothacher, J.: Stream discharge in gaged watersheds at the Andrews Experimental Forest, 1949 to present., Long-Term Ecol. Res. For. Sci. Data Bank [Database], https://doi.org/10.6073/pasta/c85f62e9070a4ebe5e455190b4879c0c, 2019.

Karlsen, R. H., Bishop, K., Grabs, T., Ottosson-Löfvenius, M., Laudon, H., and Seibert, J.: The role of landscape properties, storage and evapotranspiration on variability in streamflow recessions in a boreal catchment, J. Hydrol., 570, 315–328, https://doi.org/10.1016/j.jhydrol.2018.12.065, 2018.

Kirchner, J. W.: Catchments as simple dynamical systems: Catchment characterization, rainfall-runoff modeling, and doing hydrology backward, Water Resour. Res., 45, 1–34, https://doi.org/10.1029/2008WR006912, 2009.

Liu, J., Han, X., Chen, X., Lin, H., and Wang, A.: How well can the subsurface storage–discharge relation be interpreted and predicted using the geometric factors in headwater areas?, Hydrol. Process., 30, 4826–4840, https://doi.org/10.1002/hyp.10958, 2016.

Lyon, S. W., Koutsouris, A., Scheibler, F., Jarsjö, J., Mbanguka, R., Tumbo, M., Robert, K. K., Sharma, A. N., and van der Velde, Y.: Interpreting characteristic drainage timescale variability across Kilombero Valley, Tanzania, Hydrol. Process., 29, 1912–1924, https://doi.org/10.1002/hyp.10304, 2015.

McMillan, H., Gueguen, M., Grimon, E., Woods, R., Clark, M., and Rupp, D. E.:
Spatial variability of hydrological processes and model structure
diagnostics in a 50 km^{2} catchment, Hydrol. Process., 28, 4896–4913,
https://doi.org/10.1002/hyp.9988, 2014.

McMillan, H. K., Clark, M. P., Bowden, W. B., Duncan, M., and Woods, R. A.: Hydrological field data from a modeller's perspective: Part 1. Diagnostic tests for model structure, Hydrol. Process., 25, 511–522, https://doi.org/10.1002/hyp.7841, 2011.

Meriö, L., Ala-aho, P., Linjama, J., Hjort, J., Kløve, B. and Marttila, H.: Snow to Precipitation Ratio Controls Catchment Storage and Summer Flows in Boreal Headwater Catchments, Water Resour. Res., 55, 4096–4109, https://doi.org/10.1029/2018WR023031, 2019.

Mutzner, R., Bertuzzo, E., Tarolli, P., Weijs, S. V., Nicotina, L., Ceola, S., Tomasic, N., Rodriguez-Iturbe, I., Parlange, M. B., and Rinaldo, A.: Geomorphic signatures on Brutsaert base flow recession analysis, Water Resour. Res., 49, 5462–5472, https://doi.org/10.1002/wrcr.20417, 2013.

Ploum, S. W., Lyon, S. W., Teuling, A. J., Laudon, H., and van der Velde, Y.: Soil frost effects on streamflow recessions in a sub-arctic catchment, Hydrol. Process., 33, 1304–1316, https://doi.org/10.1002/hyp.13401, 2019.

Roques, C., Rupp, D. E., and Selker, J. S.: Improved streamflow recession parameter estimation with attention to calculation of $-\mathrm{d}Q/\mathrm{d}t$, Adv. Water Resour., 108, 29–43, https://doi.org/10.1016/j.advwatres.2017.07.013, 2017.

Rupp, D. E. and Selker, J. S.: Drainage of a horizontal Boussinesq aquifer with a power law hydraulic conductivity profile, Water Resour. Res., 41, 1–8, https://doi.org/10.1029/2005WR004241, 2005.

Rupp, D. E. and Selker, J. S.: Information, artifacts, and noise in d*Q*∕d*t* –
*Q* recession analysis, Adv. Water Resour., 29, 154–160,
https://doi.org/10.1016/j.advwatres.2005.03.019, 2006a.

Rupp, D. E. and Selker, J. S.: On the use of the Boussinesq equation for interpreting recession hydrographs from sloping aquifers, Water Resour. Res., 42, 1–15, https://doi.org/10.1029/2006WR005080, 2006b.

Rupp, D. E. and Woods, R. A.: Increased flexibility in base flow modelling using a power law transmissivity profile, Hydrol. Process., 22, 2667–2671, https://doi.org/10.1002/hyp.6863, 2008.

Rupp, D. E., Schmidt, J., Woods, R. A., and Bidwell, V. J.: Analytical assessment and parameter estimation of a low-dimensional groundwater model, J. Hydrol., 377, 143–154, https://doi.org/10.1016/j.jhydrol.2009.08.018, 2009.

Sánchez-Murillo, R., Brooks, E. S., Elliot, W. J., Gazel, E., and Boll, J.: Baseflow recession analysis in the inland Pacific Northwest of the United States, Hydrogeol. J., 23, 287–303, https://doi.org/10.1007/s10040-014-1191-4, 2015.

Santos, A. C., Portela, M. M., Rinaldo, A., and Schaefli, B.: Estimation of streamflow recession parameters: New insights from an analytic streamflow distribution model, Hydrol. Process., 33, 1595–1609, https://doi.org/10.1002/hyp.13425, 2019.

Selker, J. S. and Haith, D. A.: Development and Testing of Single-Parameter Precipitation Distributions, Water Resour. Res., 26, 2733–2740, https://doi.org/10.1029/WR026i011p02733, 1990.

Shaw, S. B. and Riha, S. J.: Examining individual recession events instead of a data cloud: Using a modified interpretation of dQ/dt-Q streamflow recession in glaciated watersheds to better inform models of low flow, J. Hydrol., 434–435, 46–54, https://doi.org/10.1016/j.jhydrol.2012.02.034, 2012.

Stewart, M. K.: Promising new baseflow separation and recession analysis methods applied to streamflow at Glendhu Catchment, New Zealand, Hydrol. Earth Syst. Sci., 19, 2587–2603, https://doi.org/10.5194/hess-19-2587-2015, 2015.

Stoelzle, M., Stahl, K., and Weiler, M.: Are streamflow recession characteristics really characteristic?, Hydrol. Earth Syst. Sci., 17, 817–828, https://doi.org/10.5194/hess-17-817-2013, 2013.

Szilagyi, J., Gribovszki, Z., and Kalicz, P.: Estimation of catchment-scale evapotranspiration from baseflow recession data: Numerical model and practical application results, J. Hydrol., 336, 206–217, https://doi.org/10.1016/j.jhydrol.2007.01.004, 2007.

Tashie, A., Scaife, C. I., and Band, L. E.: Transpiration and Subsurface Controls of Streamflow Recession Characteristics, Hydrol. Process., 33, 2561–2575, https://doi.org/10.1002/hyp.13530, 2019.

Thomas, B. F., Vogel, R. M., and Famiglietti, J. S.: Objective hydrograph baseflow recession analysis, J. Hydrol., 525, 102–112, https://doi.org/10.1016/j.jhydrol.2015.03.028, 2015.

Troch, P. A., de Troch, F. P., and Brutsaert, W.: Effective Water Table Depth to Describe Initial Conditions Prior to Storm Rainfall in Humid Regions, J. Water Resour. Res., 29, 427–434 1993.

USGS: National Water Information System, available at: https://waterdata.usgs.gov/nwis/uv?site_no=14161500 (last access: 1 February 2020), 2019.

Vannier, O., Braud, I., and Anquetin, S.: Regional estimation of catchment-scale soil properties by means of streamflow recession analysis for use in distributed hydrological models, Hydrol. Process., 28, 6276–6291, https://doi.org/10.1002/hyp.10101, 2014.

Vogel, R. M. and Kroll, C. N.: Regional Geohydrologic-Geomorphic Relationships for the Estimation of Low-Flow Statistics, Water Resour., 28, 2451–2458, 1992.

Wang, D.: On the base flow recession at the Panola Mountain Research Watershed, Georgia, United States, Water Resour. Res., 47, 1–10, https://doi.org/10.1029/2010WR009910, 2011.

Yeh, H. and Huang, C.: Evaluation of basin storage-discharge sensitivity in Taiwan using low‐flow recession analysis, Hydrol. Process., 1–14, https://doi.org/10.1002/hyp.13411, 2019.