Estimation of the base flow recession constant under human interference

The base flow recession constant, Kb, is used to characterize the interaction of groundwater and surface water systems. Estimation of Kb is critical in many studies including rainfall‐runoff modeling, estimation of low flow statistics at ungaged locations, and base flow separation methods. The performance of several estimators of Kb are compared, including several new approaches which account for the impact of human withdrawals. A traditional semilog estimation approach adapted to incorporate the influence of human withdrawals was preferred over other derivative‐based estimators. Human withdrawals are shown to have a significant impact on the estimation of base flow recessions, even when withdrawals are relatively small. Regional regression models are developed to relate seasonal estimates of Kb to physical, climatic, and anthropogenic characteristics of stream‐aquifer systems. Among the factors considered for explaining the behavior of Kb, both drainage density and human withdrawals have significant and similar explanatory power. We document the importance of incorporating human withdrawals into models of the base flow recession response of a watershed and the systemic downward bias associated with estimates of Kb obtained without consideration of human withdrawals.


Introduction
[2] Numerous studies suggest that human impacts to ecosystem degradation, climate change, increased population, and urbanization must be considered as part of the hydrologic cycle [Vorosmarty et al., 2000;Wagener et al., 2010;Vogel, 2011]. Hydrologic investigations that ignore the influence of humans have been shown to perform satisfactorily in relatively undisturbed watersheds [Mendoza et al., 2003]; however, recent research indicates that even small disturbances in hydrologic systems, such as human withdrawals [Thomas, 2012;Wang and Cai, 2009] and land-use/land cover changes [Schilling and Libra, 2003], can be important, especially in studies which focus on the low flow response of a watershed which is strongly influenced by groundwater discharge [Wang and Cai, 2010;Schilling and Zhang, 2006;Brandes et al., 2005a, among many others]. Theis [1940] argued that any newly introduced groundwater withdrawal must be balanced by one or more of the following: (1) an increase in the recharge to an aquifer, (2) a decrease in the discharge from a groundwater system to either a surface water system or a deeper groundwater system, or (3) a loss of storage to the aquifer. Given the water balance component shifts predicted by Theis [1940], it follows that increases in projected water demands require improvements in our ability to characterize low flow behavior of rivers for conjunctive water use management, maintenance of water quality, and ecosystem services. Despite the overwhelming evidence of water withdrawals and other human-induced changes impacting hydrology, most hydrologic textbooks ignore the impacts of humans on the hydrologic cycle.
[3] Hall [1968], Tallaksen [1995], and Smakhtin [2001] review the literature on the characterization of low flow behavior of a river, referred to here as base flow recession analysis. The theory of base flow recession analysis emerged from early studies of groundwater flow [Dupuit, 1863;Boussinesq, 1877;Maillet, 1905] and has since led to multiple approaches to characterize the relationship between groundwater and surface water during low flow periods. Studies which have sought to characterize coupled groundwater and surface water (gw-sw) systems normally assume a power-law relationship between base flow, Q, and groundwater storage, S. Such studies can be grouped into two categories [Dooge, 1983]: (1) linear behavior in which Q is a linear function of S, and (2) nonlinear behavior in which Q is a nonlinear function of S.
[4] Considerable controversy exists in the hydrologic literature over the issue of whether or not watersheds behave as linear or nonlinear reservoirs in the context of their low flow response. Wang and Cai [2009] illustrate the nonlinear reservoir response caused by changes in water withdrawals and the dynamics of natural watershed processes. Botter et al. [2009] derived theoretical streamflow probability density functions which occur due to the degree of linearity in storage-discharge behaviors. Watersheds that exhibit linear reservoir behavior also exhibit a fixed time constant which is often termed the base flow recession constant, which we term K b . Despite the continuing controversy over whether or not the base flow system response is linear or nonlinear [see, e.g., Wittenberg and Sivapalan, 1999], the estimation of K b continues to be required for many hydrologic models, too numerous to list exhaustively (e.g., HEC-HMS, HSPF). Moreover, the dynamics of rainfall-runoff models [Ferket et al., 2010;Beven et al., 1995] are dependent on estimation of K b or its equivalent. Vogel and Kroll [1992], Kroll et al. [2004], and Eng and Milly [2007] advocate inclusion of K b in models for estimation of low flow statistics at ungaged sites. Eckhardt [2008Eckhardt [ , 2011 illustrates that partitioning of streamflow into base flow and direct flow components generally relies on watershed storage behavior and further illustrates through sensitivity analysis how base flow separation critically depends on estimation of K b . Lo et al. [2010] provide an example of hydrologic model calibration and validation employing base flow separation to obtain runoff and groundwater discharges. Given the ubiquitous use of K b in so many hydrologic applications, it is important to (1) accurately characterize gw-sw behavior using estimates of K b and (2) characterize how the behavior of K b estimates relate to physical and anthropogenic influences known to alter hydrologic systems, the two central goals of this study.
[5] Previous investigations have explored the behavior of the base flow response as a function of watershed characteristics and climatic factors. Brandes et al. [2005b] found that drainage density and geologic indices could explain 70-80% of observed variability in the recession constant from 24 watersheds in Pennsylvania. Their results support conclusions by Zecharias and Brutsaert [1988a] who found that drainage density along with basin slope and perennial stream length explained most of the variation in base flow behavior. Zecharias and Brutsaert [1988b] noted that watershed characteristics such as basin slope, defined as the inclination angle of the impermeable layer underlying a hillslope, and drainage density explained greater variability in K b as compared to other watershed characteristics such as aquifer area percentages ( % ) and spatially heterogeneous properties such as hydraulic conductivity (k) and aquifer porosity (f). Farvolden [1963] identified a similar geologic and geomorphologic dependence to K b which was found to vary with stream channel depth and drainage density while being inversely proportional to basin relief ratios, a nondimensional ratio developed by Schumm [1954] which provides information concerning the steepness of a basin. Farvolden [1963] also identified a weak correlation between K b and the aridity ratio for arid watersheds. Likewise, van Dijk [2010] documented the humidity index as being an important factor in explaining the behavior of K b for dry catchments.
[6] Unlike previous approaches which have evaluated our ability to estimate and understand the base flow recession relationship while ignoring human withdrawals, we introduce new procedures that incorporate the impact of groundwater withdrawals on the estimation of base flow recession constants and compare their performance with traditional approaches. A primary goal of this study is to better understand the impact of human withdrawals on gwsw behavior and estimation of K b . We first consider traditional estimation methods [Barnes, 1939;Vogel and Kroll, 1992] which are given in most introductory hydrology textbooks [Dingman, 2002;Fetter, 1994;Domenico and Schwartz, 1998] and derivative approach introduced by Brutsaert and Nieber [1977]. We then consider a derivative approach from Wang and Cai [2009] which includes groundwater withdrawals. A derivation is presented for a new estimator based on the traditional approach which includes groundwater withdrawals. Finally, we employ multivariate regression procedures and the concept of elasticity to evaluate the nondimensional relative sensitivity and importance of geologic, geomorphic, climatic, and anthropogenic factors on the magnitude of K b estimates. We demonstrate that improved estimation of K b can be achieved by incorporating estimates of human withdrawals into the analysis. We also quantify the influence of human withdrawals in comparison with other watershed characteristics on observed variations in K b .

Model Development
[7] To characterize the behavior of coupled groundwater and surface water systems, we follow the work of Barnes [1939], Hall [1968], and Brutsaert and Nieber [1977], all of whom have conceptualized a simplified decay model which assumes that the initial rapid watershed input that contributed to the increase in streamflow, the rising limb of the hydrograph, have been fully depleted. It is recognized that, at sometime during the falling limb, observed streamflow is comprised entirely of discharge derived from watershed storage at which point the flow is termed base flow [Hall, 1968] and the hydrograph is termed the base flow recession curve. A base flow recession curve exhibits behavior attributed to the relationship between aquifer storage and the associated groundwater outflow to the stream channel. Hall [1968], Dooge [1983], and many others have suggested the use of a power law model to relate groundwater storage in an unconfined aquifer in direct hydraulic connection to a stream, S (mm), and the observed groundwater discharge to the stream, Q (mm d À1 ), using where and n are constants. For a linear storage-discharge relationship, n ¼ 1 in (1). Assuming a linear system, we note that in (1) has units of inverse time (t À1 ); the inverse of is often termed the time constant or recession constant [Brutsaert and Lopez, 1998;Brutsaert, 2008]. Hall [1968], Vogel and Kroll [1992], and others relate in (1) to K b where K b ¼ e À thus creating a nondimensional variable to describe base flow recession behaviors . Our attention focuses on K b instead of due to its ubiquitous use in hydrologic models. While the assumption that n ¼ 1 is not applicable for all applications, Vogel and Kroll [1992] found this to be a valid assumption in selected Northeastern United States watersheds.
[8] For an aquifer with no external transfers, continuity can be written as where I(t) represents recharge to groundwater, W(t) represents groundwater withdrawals, and ET(t) represents evapotranspiration from the groundwater table and stream. To simplify (2), we make assumptions commonly employed in previous studies. First, we assume that the analysis is conducted during stream recession periods (i.e., no precipitation events) with rapid initial percolation which ceases after precipitation ; thus, I(t) ¼ 0. Second, we assume that ET(t) has a negligible impact on (1) and (2) (1) and (2), while ignoring groundwater withdrawals, leads to where b ¼ (2n À 1)/n and a ¼ n (1/n) . Brutsaert and Nieber [1977] suggest using a log-log plot of dQ/dt versus Q to estimate the recession parameters a and b in (3). We term this approach the derivative method because Brutsaert and Nieber [1977] recommended using the time derivative of streamflow, dQ/dt, to eliminate time as the reference for estimation of base flow parameters. More recently, Wang and Cai [2009] extended (3) to include the effect of the groundwater withdrawal term in (2) resulting in [10] Combining (1), (3), and (4) leads to various estimators of K b which include the effect of groundwater withdrawals described in the next section. Though Cai [2009, 2010] were the first to introduce withdrawals into base flow analyses, this is the first study to offer several estimation methods for K b , which account for withdrawals.
2.1. Case 1: Traditional Approach (W 5 0) [11] Following Barnes [1939] and ignoring withdrawals while adhering to previous simplifying assumptions, (2) can be simplified to [12] Combining (1) and (5) while assuming Dt ¼ 1 day yields the difference equation Q tþ1 ¼ Q t K b t ; where K b is defined by Barnes [1939] as the base flow recession constant and Q t is the streamflow on day t. Defining Q 0 as the initial streamflow on day t ¼ 0, the difference equation can be linearized as where K b1 is the estimate of the base flow recession constant for Case 1. Equation (6) can be solved graphically or mathematically by linear regression approaches where ln(Q t ) is the dependent variable and t is the independent variable with intercept ln(Q 0 ) and slope equal to ln (K b1 ). An identical model assuming independent model errors with zero mean and constant variance was investigated by Vogel and Kroll [1996] using ordinary least squares regression to estimate K b1 . They showed that the use of log-space error in (6) led to improved performance of the resulting estimators compared to the use of real space errors.
2.2. Case 2: Traditional Approach (W 6 ¼ 0) [13] Here we introduce a new estimator of K b that includes the impact of human withdrawals. Including groundwater withdrawals which deplete aquifer storage, the continuity equation becomes [14] Combining (1) and (7) with k ¼ 1/ and using subscripts to represent time dependence, then [15] Here we assume groundwater withdrawals are constant during each individual recession hydrograph and are therefore not a function of time during an individual hydrograph recession. This assumption is necessary to incorporate monthly groundwater withdrawals into the analysis. Multiplication of (8) by the integration factor e t/k leads to [16] Using the chain rule and integration leads to where the integration constant is kQ 0 and K b2 ¼ e À1/k or e À . Equation (10) can be linearized by noting that [17] Equation (11) is a linear model with dependent variable ln(Q t þ W), independent variable time t, intercept ln(Q 0 þ W), and slope ln(K b2 ). Thus, ordinary least squares regression may be used to estimate the base flow recession constant in (11).

Case 3: Derivative Approach (W 5 0)
[18] Given (3) from Brutsaert and Nieber [1977] and Vogel and Kroll [1996] illustrate that the intercept parameter a is equal to ln(K b ). Taking the limit of (3) as b approaches 1 and taking logs yields [19] Equation (12) is a linear model with a slope term constrained to 1 on a plot of -dQ/dt versus Q. The estimate of K b3 is obtained by fitting (12) to individual recession data; this approach deviates from previous studies who suggest that individual recessions exhibit a slope of 1.5 or 3 [Brutsaert and Nieber, 1977] or 2 [Biswal and Marani, 2010;Shaw and Riha, 2012]. Consistent with Case 3, we assume that the system is linear and hence exhibits a slope of 1. A simple and attractive estimator of the intercept [ln(Àln(K b3 )] is the mean of the difference [ln(ÀdQ/ dt) À ln(Q t )].
[20] Brutsaert and Nieber [1977] estimated dQ/dt using the finite difference scheme given by with mean discharge Q estimated using [Q i þ Q iÀ1 ]/2; thus errors arise in both estimation of dQ/dt and the mean value of Q. Rupp and Selker [2006] summarized the effect of the choice of Dt and precision of flow estimates on recession plot analysis. Thomas [2012] recommended more robust generalized cross-validation (GCV) smoothing spline alternatives over the arbitrary numerical difference schemes for the estimation of dQ/dt. Here estimates of dQ/dt in (12) were obtained using optimal smoothing splines which were fit to each individual recession hydrograph within the Rproject [R Development Core Team, 2012]. Employing GCV methods with smoothing splines provide nearly optimal estimates of dQ/dt [Craven and Wahba, 1979] in addition to preserving the observed flow Q in contrast to methods used in most previous recession studies.
2.4. Case 4: Derivative Approach (W 6 ¼ 0) [21] Taking logarithms to linearize (4) and applying the limit as b approaches 1 yields [22] As done for Case 3, the time derivative dQ/dt was estimated using smoothing splines fit to individual hydrograph recessions. An estimate of K b4 can be obtained from equation (13) by noting that the mean of the difference [ln(ÀdQ/dt) À ln(Q t þ W)] is equal to the mean value of ln(Àln[K b4 ]), similar to the estimation scheme for Case 3.

Database Description
[23] We apply our base flow estimation schemes to watersheds in New Jersey, USA, because it is currently one of the only states in which time series of monthly groundwater withdrawals are available. A total of 15 watersheds were selected for this study because (1) these watersheds are within high density population areas of northern New Jersey which meets the challenge of Sivapalan et al. [2011] and Vogel [2011] to study anthropogenic impacts on hydro-logic processes where humans live and (2) a long-term ($19 year) monthly database of groundwater withdrawals is publicly available [New Jersey Geological and Water Survey (NJGWS), 2011]. Daily streamflow records at U.S. Geological Survey (USGS) stream gages were used to isolate suitable recession hydrographs for analysis ( Figure 1) using the hydrograph selection algorithm described by Vogel and Kroll [1996]. That algorithm identifies the beginning of a streamflow recession when a 3 day moving average begins to decrease and ends when a 3 day moving average begins to increase. Recessions with lengths greater than or equal to 10 days were used for this study. Additionally, the first three points of the recession were removed for estimation of the base flow recession constant to remove the influence of other runoff processes. Vogel and Kroll [1992], Brutsaert and Lopez [1998], and others have identified different recession hydrograph isolation algorithms ; Stoelzle et al. [2013] identified the approach of Vogel and Kroll [1992] as accurately characterizing average recession behaviors compared to algorithms used by Brutsaert and Nieber [1977] and Kirchner [2009], though their study raises questions regarding the uniqueness of any base flow recession characterization. USGS stream gage sites were selected based on the continuous daily records for the calendar years 1990-2009 to coincide with available withdrawal data. Available geographical information databases obtained from the New Jersey Department of Environmental Protection Bureau of Geographic Information Systems (GIS) were used to select watersheds ; for the purpose of this study, we selected watersheds with georeferenced surficial aquifers which did not intersect multiple watershed boundaries and which exhibited possible hydraulic connection to streams based on the visual intersection of aquifer and hydrography layers (Figure 1).
[24] Brutsaert [2005] and Vogel and Kroll [1992] derive a relationship between the base flow recession constant and aquifer properties; thus, the geologic and hydrogeologic characterization of the region theoretically dictates values of K b . Watersheds selected for this study are located in northern New Jersey within the Valley and Ridge, Highlands, and Piedmont geologic provinces characterized as having valley fill aquifers [Kummel, 1940]. The Valley and Ridge Province is underlain by faulted and folded sedimentary bedrock ; such structural characteristics have undergone erosion resulting in removal of sandstone creating valleys and resistant shale/limestone ridges. The Highlands region is underlain by igneous and metamorphic bedrock ; these bedrock types are resistant to weathering resulting in deeply incised and steep stream valleys. The Piedmont Province is underlain by both igneous and sedimentary bedrock creating a broad lowland area. The unconsolidated formations that comprise aquifers in northern New Jersey are characterized as glacial Pleistocene stratified drift from the Wisconsin ice sheet or from more recent valley fill deposits [Kummel, 1940]. A summary of physical watershed characteristics for the selected watersheds illustrated in Figure 1 is included in Table 1; a full description of data sources for this study is included in section 5.1.
[25] Daily precipitation data were obtained from the National Climatic Data Center (NCDC) and the Utah Climate Center database for sites with nearly continuous daily precipitation records for the period 1990-2009 while monthly data were obtained from PRISM [Daly et al., 1994]. Selected daily precipitation sites are illustrated in Figure 1. For the purpose of this study, precipitation totals reported as trace (T) are considered to be zero. Precipitation data were used to eliminate any recession hydrograph during which streamflow was observed to decrease despite a concurrent precipitation event totaling greater than 0.10 mm measured at the two precipitation gage sites nearest to the streamflow gage.
[26] Deterministic approaches for the determination of the impact of well withdrawals on surface water flows [e.g., Hunt, 2012] illustrate spatially dependent relationships  where large withdrawals in close proximity to a stream result in a larger impact on streamflow as compared to a small withdrawals located a great distance from a stream. Further, withdrawals from unconfined aquifers with a potential hydraulic connection to the stream will produce larger impacts than withdrawals from confined aquifers where the confining unit is located between the screened material and the stream. Given these concerns, only groundwater withdrawals from unconfined aquifers were used from available time series of monthly groundwater withdrawals over the period 1990-2009 obtained from the New Jersey Geological Survey [NJGWS, 2011]. Total watershed withdrawals were calculated by aggregating HUC14 withdrawals within delineated watersheds. Total watershed withdrawals were assumed to be uniformly distributed across each month to obtain daily withdrawal rates.

Experimental Design
[27] The four K b estimators summarized in Table 2 are applied at each site i resulting in m i estimates of K b where m i is the number of hydrograph recessions at site i. Seasonal variability in gw-sw behaviors has been noted in previous studies [Szilagyi et al., 2007;Kroll, 1989]. To characterize seasonal effects on the estimation of K b , recession data from each site was separated into four seasons, winter (December-February), spring (March-May), summer (July-August), and fall (September-November), where each season was assigned based on the last day of the hydrograph recession.

Results: K b Estimation
[28] Vogel and Kroll [1996] observed less variability in estimates of the base flow recession constant for longer recessions. To characterize each site and season with a single average estimate of K b , a recession-length weighted average was employed : where t i,s,j is the total number of days of the recession hydrograph period used to estimate the base flow recession constant, j is the number of recession hydrographs and K bi,s,j is the estimate of K b from the jth hydrograph recession at site i during season s. Figure 2 illustrates recessionlength weighted averages of K b from (14) by both season and estimation scheme. The results in Figure 2 illustrate that the lowest and most variable estimates are obtained from Case 1 as compared to other estimation schemes (where number of estimates, n, is 59, 59, 58, and 58, respectively, for Cases 1, 2, 3, and 4). Seasonal estimates obtained from Cases 3 and 4, the derivative approaches, illustrate differences in estimates of K b in low ET seasons (spring and winter) as compared to high ET seasons (summer and fall). Previous studies documented that much of the variability in base flow recessions between seasons depended on evapotranspiration [Szilagyi, 1999;Wittenberg, 2003] or the active drainage network [Biswal and Marani, 2010]. Figure 2 illustrates that greater variance in estimators was identified during seasons which exhibit higher ET. For all estimators tested, estimates of K b obtained from summer recessions exhibited the lowest variance compared to other seasons. For all seasons, K b estimates obtained from Cases 3 and 4 exhibited lower variance when compared to estimates obtained from Cases 1 and 2. Such results agree with Szilagyi et al. [2007] and Kirchner [2009] who suggest that reduced ET should reduce variability in observed base flow recession parameters. Additional variability in base flow recession constant estimation may be due to varying hydrogeologic properties of the aquifers within a watershed across different flow regimes as described by previous studies [Brutsaert and Nieber, 1977;Brutsaert and Lopez, 1998;Brutsaert, 1988a, 1988b;Farvolden, 1963]. Further, groundwater elevations tend to be highest in the spring and lowest in the fall as evidenced in monthly groundwater elevation data obtained from USGS across New Jersey. Such changes in groundwater elevations may alter active stream networks as described by Biswal and Marani [2010]. As described by Van de Griend et al. [2002], physical constraints of aquifer systems such as the decrease in hydraulic conductivity with depth in stratified drift aquifers due to compaction can greatly affect the seasonal storage-discharge behaviors.
[29] Here we introduced a dimensionless withdrawal variable to describe the degree of withdrawal. The variable, , characterizes the impact of groundwater withdrawals as a ratio to overall streamflow during the recession where w j is the daily average groundwater withdrawal for hydrograph recession j, Q r j represents the average streamflow during the jth hydrograph recession and again m i is the total number of observed hydrograph recessions available at site i.
[30] Figure 3a illustrates a comparison of estimates of K b from individual hydrograph recessions obtained from equation (6) (Case 1) and equation (11) (Case 2) for different ranges of . Figure 3a (n ¼ 827) illustrates that some variability in estimates of K b obtained from individual hydrograph recessions can be attributed to the impact of groundwater withdrawals; this is most clearly illustrated by estimates during periods in which is high (points representing 0.10<<0.25) which show a range in estimated K b values from 0.80 to 0.97 obtained for Case 1 with a more limited range from 0.85 to 0.97 for Case 2. It must be stressed, however, that K b is bounded (K b 1) and thus larger estimates of K b will naturally result in lower variability. Figure 3b (n ¼ 817) compares estimates of K b for Case 3 versus Case 4. Figure 3b illustrates a similar systematic impact of groundwater withdrawals on estimates of the base flow recession constant. Figure 3 illustrates that incorporation of groundwater withdrawals into estimation of the base flow recession constant generally results in an increased value of K b as evidenced the fact that K b2 >K b1 ( Figure 3a) and K b4 >K b3 (Figure 3b). A significant   (14) for cases summarized in Table 2. Statistical differences between estimates were tested using the Wilcoxon signed rank sum test [Helsel and Hirsch, 2002]. Statistical differences with ¼ 0.05 were identified between summer/fall and winter/spring in (c) and (d), between Cases 1 and 2 for all seasons (e-h), and between Cases 1 and 3 and between Cases 1 and 4 during summer and fall (g and h). Boxplots represent 75% and 25% quantiles while whiskers represent 90% and 10% declies. difference in estimates of K b1 and K b2 can be inferred from standard error estimates of the base flow recession constants for Cases 1 and 2. Overall it was found that withdrawals exceeding, on average, values of ¼ 0.084 resulted in significantly different estimates of K b1 and K b2 . We conclude that when an average of 8.4% or more of streamflow is withdrawn during a recession event such withdrawals will significantly impact our ability to characterize recession behaviors without considering human interferences.

Relating Base Flow Recession Constants to Watershed Characteristics
[31] In this section, we attempt to describe the relationship between K b and physical watershed characteristics and human interferences to characterize the relation between base flow recession and its controlling parameters. Multivariate regression has been used to relate the base flow recession constant to physical [Brandes et al., 2005b;Brutsaert, 1988a, 1988b] and climatological [Farvolden, 1963;van Dijk, 2010] variables. Our hypothesis is that K b is a function of physical, anthropogenic, and climatic influences and, therefore, is no longer dependent solely on physical characteristics of the watershed. Our approach to evaluate this hypothesis includes development of a multiple regression model that accounts for the primary physical determinants related to K b as well as those due to climatic and anthropogenic influences.
[32] Vogel and Kroll [1992] provide a physically based model relating K b to watershed characteristics using simplifying watershed and hydrogeologic assumptions (including Dupuit flow) as where k is the average watershed hydraulic conductivity, H is the basin relief, d d is the drainage density, % is the fraction of the watershed underlain by aquifers contributing to base flow and f is the drainable soil porosity, sometimes referred to as the effective porosity or specific yield.
[33] The physical model in (16) relates the base flow recession constant to physical watershed characteristics. Our results suggest that K b is not a constant related solely to physical characteristics of a watershed and, therefore, may be dependent on other exogenous variables including seasonality and groundwater withdrawals.

Multivariate Regression Database
[34] Watershed characteristics needed for the proposed multivariate models include, at a minimum, aquifer permeability, basin elevations, porosity, the proportion of the watershed underlain by surficial aquifers, and the drainage density. Additional factors including hydroclimatic indices, a seasonality index, and the groundwater withdrawal index from (15) were also included in model testing. Availability of hydrogeologic information including hydraulic conductivity and porosity are limited ; thus, lumped estimates of porosity were assigned based on the reported sedimentary materials while transmissivity, a parameter related to hydraulic conductivity and aquifer depth, was used given the availability of groundwater pumping test data. Limitations of data used in the analysis are included below. A summary of potential explanatory variables is shown in Table 3.

Transmissivity (T)
[35] Most watersheds used as part of this study have reported estimates of transmissivity derived from groundwater pumping tests [New Jersey Department of Environmental Protection (NJDEP), 2012]. Estimates of transmissivity are point estimates; therefore, area-weighted averages were used via Theissen polygons to obtain representative watershed transmissivity. We recognize that such methods do not provide representative estimates of bulk watershed properties; however, for this initial analysis to characterize the effect of withdrawals on gw-sw behaviors, we used available data which included point estimates and aquifer characteristics. For watersheds without estimates of transmissivity, methods were used to relate watershed properties with transmissivity estimates to those without documented transmissivity; such sites are indicated in Table 1. In general, aquifers with high aquifer productivity and characterized as being composed of sand and gravel exhibited high transmissivity estimates (T > 500 m 2 /d) while lower aquifer productivity and characterization as till or sand/silt exhibited low transmissivity (T < 300 m 2 /d).

Alpha (a % )
[36] The aerial percentage of the watershed underlain by aquifers was calculated using GIS data available through NJGIS [Herman et al., 1998]. Estimates of % represent the ratio of area of aquifers divided by the total contributing watershed area.

Porosity (f)
[37] Average watershed porosity was assigned based on the available GIS coverages [Herman et al., 1998] describing aquifer productivity and aquifer formations and linked to porosity ranges (0.1, 0.2, and 0.3) obtained from Fetter [1994]. For example, an aquifer characterized as being composed of silt, sand, and clay was assigned a porosity of 0.10 while sand and gravel aquifers with high yields were assigned porosity values of 0.30. Effective porosity estimates used in this study are higher than reported estimates obtained from (16) and similar equations introduced by Brutsaert and Nieber [1977] by orders of magnitude [Brutsaert and Nieber, 1977, f ¼ 0.045;Eng and Brutsaert, 1999, f ¼ 0.017;Brutsaert and Lopez, 1998, f ¼ 0.0167;Rupp et al., 2009, specific yield $0.0021;Brutsaert, 2008,  specific yield $0.024]. The potential ranges of effective porosity have been described by Johnson [1968]. Consider an aquifer unit is described as silty-sand/clay sand; Johnson [1968] illustrates that specific yield can vary from 0.04 to 0.30. Further, complicating estimates used in this study is that production well testing is likely to occur in the most favorable regions of aquifers and will not capture heterogeneities in aquifer systems known to exist.

Drainage Density (d d )
[38] The New Jersey Department of Environmental Protection (NJDEP) provides National Hydrography Dataset (NHD) GIS data which includes NHDFlowline feature classes, or ''blue lines'' representing streams. The GIS database was used to compile watershed variables needed to calculate drainage density which we define as the ratio of the total stream channel length to drainage area following Vogel and Kroll [1992].

Degree of Withdrawal (c)
[39] The variable , defined in (15), was used to characterize the degree of groundwater withdrawals which are expected to impact estimation of K b . 5.1.6. Average Watershed Slope (S) [40] Previous studies employed the concept of basin relief (H) which Zecharias and Brutsaert [1985] define as the difference between the basin summit elevation and the basin outlet elevation. Vogel and Kroll [1992], for example, estimate average watershed slope, S, using the relation S % 2Hd d based on the work by Strahler [1950]. For this study, average watershed slope was calculated using the approach by the Natural Resources Conservation Service (NRCS) which calculates average watershed slope by [41] Contour lines were created in GIS from New Jersey digital elevation grid data [New Jersey Geological Survey (NJGS), 1999]. 5.1.7. Runoff Ratio (RR) [42] The seasonal runoff ratio, the ratio of seasonal streamflow to seasonal precipitation (Q/P) defined by Chang [2007], characterizes the seasonal behavior of watershed response to precipitation. The runoff ratio (RR) is a function of climate and watershed properties as it depends on antecedent moisture conditions and maximum soil water retention. For this study, monthly precipitation data were obtained from PRISM and area-weighted averages for each watershed. PRISM data were selected for this study since the spatial regression procedure to produce the spatial and temporal data set accounts for topographic impacts on precipitation [Daly et al., 1994]. Average monthly streamflow was obtained from USGS stream gages and normalized by contributing watershed area. 5.1.8. Aridity Index (AI) [43] The seasonal aridity index (P/PET) defined as the ratio of seasonal precipitation (P) to seasonal potential evapotranspiration (PET) compares available water to the energy within a system. To calculate the seasonal aridity index, monthly precipitation data was obtained using watershed-area-weighted averages from PRISM [Daly et al., 1994] while PET was estimated using the Hargreaves method [Shuttleworth, 1993].

Palmer Drought Stress Index (PDSI)
[44] PDSI was introduced by Palmer [1965] to evaluate soil moisture deficiency by accounting for the interaction of precipitation and evapotranspiration. Alley [1984] documented deficiencies with PDSI including its inability to characterize temporal estimation of water budget components. In spite of this criticism, PDSI remains commonly used as a climate index [e.g., Huang et al., 2011]. Monthly PDSI was obtained for the regions of New Jersey from the National Climatic Data Center (NCDC).

Seasonality
[45] Estimates of K b are shown to vary depending on season by Kroll [1989], Szilagyi et al. [2007], and Aksoy et al. [2001] and as illustrated in Figure 2. To account for the potential relationship between K b and season, a seasonality term is introduced as where ! is (2/4) and t is the season numbered from 1 to 4 (winter-fall) following Helsel and Hirsch [2002].

Model Development
[46] Ordinary least squares (OLS) multivariate regression procedures were used to estimate the model parameters of the hybrid physical/statistical model of the form given in (16). To enable the use of multivariate OLS regression we fit the model where X i , . . . , X n are the natural logs of the explanatory variables summarized in Table 3, 0 , . . . , n are model coefficients and " i are normally distributed model errors with zero mean and constant variance 2 . As summarized in Table 3, potential explanatory variables include climatic and anthropogenic variables in addition to physical watershed attributes. The dependent variable for (19) is obtained from recession-weighted averages from (14) and thus represents seasonal estimates of K b from all sites for each case. Various multivariate regression model selection methods were used including stepwise regression analysis and best subsets regression in addition to the working knowledge of the physical model form for K b in (16). For this study, we only include independent variables whose model coefficients are significantly different from zero using a 5% level t test. In addition, influence statistics were calculated to identify observations that exhibited unrealistic influence on the regression model parameter estimates. To ensure that statistical inference associated with our model testing was reasonable, we performed diagnostics to ensure the normality, independence, and homoscedasticity of the model residuals. Finally, variance inflation factors were examined to assure against multicollinearity among the explanatory variables. In general, we followed the guidelines of Helsel and Hirsch [2002] for the diagnostic evaluation of the resulting multivariate models summarized in Table 4.
[47] For all cases, significant explanatory variables summarized in Table 4 include drainage density (d d ), average watershed slope (S), transmissivity (T), groundwater withdrawals (), and aquifer percentage ( % ). One would expect the variable season to be significant given results in Figure 2; however, one must note that is a function of streamflow and withdrawals (15) which are known to depend on season. Thus, a seasonal signal is already inherent in the model via incorporation of the impact of water withdrawals. Climate variables were not found to be significant for the sites considered; such a result may be due to the limited spatial area employed for this study over which many climatic variables are relatively homogeneous.

Evaluation of Multivariate Model Results
[48] McCuen et al. [1990] and many others have documented retransformation bias due to retransformation of a model which was fit using a log transformation of model variables, as is the case in (19). The bias was calculated for each model in Table 4; the relative bias was found to be low (maximum relative bias of 0.009 for the four cases), thus no bias correction factors were employed. Nash-Sutcliffe efficiency [Nash and Sutcliffe, 1970] was calculated for estimated values of K b and is presented in Figure 4. Figure 4 illustrates a comparison between estimated K b from the multivariate regression based on the recession-weighted averages at each site and season. We avoid comparisons of R 2 between models due to the potential for spurious correlation for Cases 2 and 4 which included in the estimation scheme and in the multivariate regression results. Results summarized in Table 4 and Figure 4 indicate that the multivariate models which include groundwater withdrawals exhibit a lower prediction sum of squares (PRESS). The PRESS statistic is defined by P nÀ1 i¼1 e 2 i where e i is the prediction (or delete 1) residual (see Helsel and Hirsch [2002] for further background on prediction residuals). Here PRESS units are the square of the dependent variable which in our case is dimensionless. Further, we note that Case 2 exhibits the lowest PRESS and highest NSE among the cases considered. We conclude, for this particular application, that the traditional semilog base flow recession approach adapted to include groundwater withdrawals yields more physically realistic estimates of K b as compared to the other estimators tested.
[49] To evaluate the sensitivity of estimates of K b to changes in the watershed variables used to explain their variability, we employ the nondimensional sensitivity concept of elasticity. The units of explanatory variables (Table  4) are not identical ; thus, it is difficult to interpret their resulting impact on K b . Elasticity is an attractive concept because it is a dimensionless approach to evaluating sensitivity defined as where x is an explanatory variable used to predict K b as summarized in Table 4. Elasticity " x represents the percentage increase (decrease) in the magnitude of K b , which results from a one percent increase (decrease) in the explanatory variable x. Here @K b @x is estimated from the regression model in (19) for each case and evaluated for each site and for each season. We define elasticity of: geomorphology (d d ), geology (T and % ), topography (S), and human withdrawals () of K b . The absolute value of the elasticity of each explanatory variable is shown using boxplots in Figure 5 to illustrate the relative importance of each variable. For all cases illustrated in Figure 5, the geomorphology elasticity of the base flow recession constant was found to be the most influential factor impacting K b . This result is consistent with previous research which suggests that drainage density is the main contributor to spatial variations in K b [Zecharias and Brutsaert, 1988b, Brandes et al., 2005b, Farvolden, 1963Brutsaert and Nieber, 1977]. Importantly, the next most influential factor for all cases was the human withdrawal elasticity of K b . For Cases 1 and 2, the value of the elasticity of human withdrawals is roughly equal to the geomorphology elasticity. This result highlights the high sensitivity of K b to variations in human groundwater withdrawals.
[50] Figure 5 shows that the elasticity of the remaining explanatory variables was small, especially when compared to the water withdrawal and geomorphology elasticity. These results support the results of previous research which documented the influence of basin slope (S) Brutsaert, 1988a, 1988b;Farvolden, 1963], aquifer area ( % ) [Brutsaert and Nieber, 1977], and aquifer conductivity (T) [Brutsaert and Nieber, 1977] on base flow recessions. However, it must be stressed that K b itself is a sensitive parameter ; a change in K b from 0.93 to 0.92, for example, is a large difference and can greatly affect results such as low flow estimation [Kroll et al., 2004] and base flow separation procedures [Eckhardt, 2011].

Conclusions
[51] This study has explored the sensitivity of the behavior of groundwater/surface water interactions using the base flow recession constant, K b , which has widespread applications in hydrologic engineering and science. We have summarized the application of K b in testing hydrologic science hypotheses and its widespread use in hydrologic models. Given the complex coupling of human and natural hydrologic processes both temporally and spatially, there is a great need to incorporate anthropogenic influences into hydrologic models in general, and in the estimation of base flow recession constants in particular. This study sought to accurately characterize groundwater-surface  (19)) along with their p values (in parentheses). All coefficients were significant at the 5% level.
water (gw-sw) behavior and characterize how the behavior of K b relates to both physical and anthropogenic influences. We have shown that when groundwater withdrawals are ignored, characterization of gw-sw behavior can be systematically biased, thus confounding our ability to model such systems. As increased attention is devoted to groundwater storage to meet water demands, this study illustrates the necessity to incorporate anthropogenic effects, in this case groundwater withdrawals, to properly characterize base flow recessions.
[52] One of the primary goals of this study was to develop a new approach for integrating human withdrawals into the estimation of the base flow recession constant. Though Wang andCai [2009, 2010] were the first to introduce withdrawals into base flow analyses, this is the first study to introduce and compare the performance of several estimation methods for K b , which account for water withdrawals. Our approach extends traditional semilog methods and derivative methods to incorporate the influence of human withdrawals. Our findings indicate that approaches which include human withdrawals to obtain the base flow recession constant (Cases 2 and 4) result in systemically higher estimates of K b than traditional approaches which ignore groundwater withdrawals. Furthermore, for all cases considered in this study, the traditional semilog method for estimating K b modified to include human withdrawals (Case 2) exhibited seasonal variability ( Figure 3) and led to multivariate models for predicting K b with higher precision (Table 4) than all other estimators considered. Our findings indicate that traditional approaches adapted to account for water withdrawals yield more consistent (Figure 3) and physically realistic (Figure 4 and Table 4) estimates of K b than traditional approaches which ignore withdrawals as well as methods which employ derivative approaches. Our results illustrate that average withdrawals of roughly 8% or more of streamflow during a recession event impact our ability to characterize recession behaviors.
[53] We also sought to quantify the role of geologic, geomorphic, topographic, and anthropogenic factors on variations in the value of K b . For this purpose, the concept of dimensionless elasticity is used. Geomorphologic (drainage density) and anthropogenic (human withdrawal) elasticities were found to be the most influential factors impacting the base flow recession constant, K b . Previous studies which ignored anthropogenic influences (human withdrawals) suggested that geomorphic variables such as drainage density or topographic variables such as basin slope are the most important variables for explaining variations in K b ; our study supports those previous findings but illustrates that groundwater withdrawals can exhibit an equally significant and important influence on observed variations in K b .  Table 4. Boxplots represent 75% and 25% quantiles while whiskers represent 90% and 10% declies.