Assessing the resiliency of surface water and groundwater systems under groundwater pumping

Supplement for Assessing the resiliency of surface water and groundwater systems under groundwater pumping S. B. Seo1, G. Mahinthakumar2, A. Sankarasubramanian2 and M. Kumar3 1Department of Civil and Environmental Engineering, Institute of Engineering Research, Seoul National University, Seoul, Republic of Korea. 2Department of Civil, Construction, and Environmental Engineering, North Carolina State University, Raleigh, North Carolina, USA. 3Nicholas School of Environment, Duke University, Durham, North Carolina, USA.


Abstract:
Since surface water and groundwater systems are fully coupled and integrated systems, increased groundwater withdrawal during drought may reduce groundwater discharges into the stream, thereby prolonging both systems' recovery from drought. To analyze watershed response to basin-level groundwater pumping, we propose an uncertainty framework to understand the 5 resiliency of groundwater and surface water systems using a fully-coupled hydrologic model under transient pumping. The proposed framework incorporates uncertainties in initial conditions to develop robust estimates of restoration times of both surface water and groundwater systems and quantifies how pumping impacts state variables such as soil moisture. Groundwater pumping impacts over a watershed were also analyzed under different pumping volumes and different 10 potential climate scenarios. Our analyses show that groundwater restoration time is more sensitive to variability in climate forcings as opposed to changes in pumping volumes. After the cessation of pumping, streamflow recovers quickly in comparison to groundwater, which has higher persistence. Pumping impacts on various hydrologic variables were also discussed. Given that surface water and groundwater are inter-connected, optimal management of the both 15 resources should be considered to improve the watershed resiliency under drought. Potential for developing optimal conjunctive management plans using seasonal-to-interannual climate forecasts is also discussed.

Introduction
Groundwater is an important source of water for various needs, including public supply, agriculture, and industry (Barlow and Leake, 2012;Sankarasubramanian et al., 2017). As the demand for groundwater use continues to rise with increasing population, groundwater resources are depleting faster than they can be replenished (Gleeson et al., 2012). Unless sustainable 5 withdrawal strategies are implemented, continued appropriation of groundwater resources can lead to streamflow depletion since groundwater and surface water are a tightly coupled interconnected system (Muller and Male, 1993;Winter et al., 1998;Sophocleous, 2002;Kumar et al., 2009b). Therefore, understanding of groundwater pumping impacts on water-budget components is critical for conjunctive management of both surface water and groundwater 10 resources.
Many studies have evaluated groundwater pumping impacts on streamflow based on conceptualized interaction between surface water and groundwater (Muller and Male, 1993;Konikow and Leake, 2014). Rasmussen et al. (2003) developed a semi-analytic integral solution for the transient response of what associated with oscillatory pumping. Recent advances in 15 computational power has also facilitated intensive simulations using extensive data sets (Kendy and Bredehoeft, 2006;Filimonova and Shtengelov, 2013;Konikow and Leake, 2014).
Forementioned studies demonstrated the dependency of streamflow on aquifer properties and distance of the pumping wells in hypothetical cases without considering the complexity of hydrogeological processes in actual watersheds. On the other hand, case studies evaluating the 20 impacts of groundwater pumping on surface water resources in real hydrologic watershed systems have generally been performed using groundwater models without being fully coupled.
For example, MODFLOW (Harbaugh, 2005) has been popularly used to simulate impacts on surface water resources due to groundwater pumping (Scibek et al., 2007;Zume and Tarhule, 2007;Garner et al., 2013). Studies also coupled a hydrologic model with a groundwater model to evaluate groundwater pumping feedbacks on surface water (Kollet and Maxwell, 2008;Lin et al., 2013;Condon and Maxwell, 2014). However, these studies examined groundwater pumping 5 impacts on future water availability at the regional scale by assuming hypothetical "periodic" pumping and sparse representation of heterogeneities in aquifer and geology structure.
Regardless of the type of simulation model (i.e., a groundwater model alone, or a coupled surface water and groundwater modeling approach), studies (Kendy and Bredehoeft, 2006;Zume and Tarhule, 2007;Barlow and Leake, 2012;Garner et al., 2013) primarily focused on estimating 10 the equilibrium state -a steady state having reduced streamflow and groundwater storageinduced by "continuous/periodic" pumping series. Further, transient groundwater withdrawal under a drought can cause adverse impacts such as depletions in streamflow and groundwater storage (Barlow and Leake, 2012;Garner et al., 2013). Studies have also assessed groundwater pumping impacts on water-budget components at the watershed scale under potential climate 15 change (e.g., Woolfenden and Nishikawa, 2014). However, these studies have mostly been carried out in dry climate regions such as Arizona and California. (e.g., Leake and Pool, 2010;Garner et al., 2013;Woolfenden and Nishikawa, 2014). Here, we focus on a humid basin from a rainfall-runoff regime in which analyze the groundwater pumping impacts conditioned on drought conditions. 20 The main objective of this study is to analyze transient pumping impacts on streamflow, groundwater and other water budgets under observed and potential climatic conditions. To account for feedback from groundwater pumping on all the hydrologic variables, we employed a Hydrol. Earth Syst. Sci. Discuss., https://doi.org /10.5194/hess-2017-402 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 10 August 2017 c Author(s) 2017. CC BY 4.0 License. fully coupled hydrologic model, Penn-state Integrated Hydrologic Model (PIHM), which simulates both surface and sub-surface hydrologic processes over prescribed triangular elements.
We primarily focused on the resiliency -the combined time needed for surface water or groundwater system to reach equilibrium state after transient pumping and then to recover to prepumping conditions after the cessation of groundwater pumping -of the watershed and also the 5 associated streamflow and groundwater depletions due to the transient pumping. Understanding and quantifying transient pumping impacts at the river basin scale could be addressed by analyzing the time series of water-budget components (e.g., soil moisture, streamflow, groundwater) obtained under "pumping" and "no-pumping" conditions. We first show how the restoration times of surface water and groundwater can vary by simple comparison of the two 10 (i.e. "pumping" and "no-pumping") time series. Since estimation of restoration times and depletion volumes could be estimated only through modeling efforts, we argue that it is important to consider relevant uncertainties. Hence, we proposed the uncertainty framework that rigorously incorporates uncertainties in initial conditions and input variables (Li et al., 2009;Sinha and Sankarasubramanian, 2013;Yossef et al., 2013). Based on the proposed uncertainty 15 framework, we developed the "null distribution of simulated streamflow" under "no-pumping" and then compare the estimated streamflow under "pumping" to understand the pumping impacts on watershed variables. We also considered the pumping impacts on watershed resiliency under different demand scenarios and climatic conditions. The proposed uncertainty framework is demonstrated for the Cape Fear watershed in North Carolina (NC) which is experiencing severe 20 shortages due to significant increase in demand over the past three decades (Singh et al., 2014).
Findings from the study are also discussed and generalized on how conjunctive management plans need to be developed considering both changing climate and demand patterns.  Figure 1). The headwaters of Haw River run 177 km into the Jordan Lake reservoir and the drainage area of the basin is 3,945 km 2 . The Haw river primarily provides freshwater for residential (cities such as Greensboro, Burlington, and Durham), industrial and recreational uses. The climate of the Haw River basin is characterized by humid subtropical climate receiving 1,138 mm of mean annual precipitation. Figure 1a presents land cover 10 classification of the Haw River basin. More than half of the entire watershed is still covered by forest.
The land surface in the piedmont region is underlain by clay-rich, unconsolidated material derived from in situ weathering of the underlying bedrock. This material, which averages about 10 to 20 m in thickness, is referred to as saprolite (Heath, 1984). Aquifers of the 15 Piedmont region are localized, complex fractured metamorphic, igneous, and sedimentary rocks.
The rocks are covered almost everywhere by regolith, and most of groundwater is stored in the shallow, porous regolith (Lindsey et al., 2006). Since unconfined groundwater generally occurs both in the porous and shallow regolith (Heath, 1984), it is assumed that groundwater is stored in unconfined aquifer across the entire watershed. The regolith contains water in pore spaces 20 between rock particles. The bedrock, on the other hand, does not have any significant intergranular porosity. The hydraulic conductivities of the regolith and the bedrock are similar and Hydrol. Earth Syst. Sci. Discuss., https://doi.org /10.5194/hess-2017-402 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 10 August 2017 c Author(s) 2017. CC BY 4.0 License. range from about 0.001 to 1 m/day (Heath, 1984). With respect to recharge conditions, it is noted that forested areas of the piedmont region have thick and very permeable soils overlain by a thick layer of forest litter. As water needs in the Piedmont region increase and as commissioning of new reservoirs is becoming increasingly difficult due to environmental impacts, it will be necessary to make more intensive use of groundwater (Heath, 1984). 5 Figure 1e shows the potential pumping locations. Multiple locations of pumping wells were hypothetically assigned for withdrawing water for public supply and irrigation purpose. A total of 150 wells for public supply were uniformly distributed in urban areas (i.e., City of Greensboro, Burlington, and Durham) and a total of 150 wells for irrigation were randomly distributed in the agricultural land use areas (i.e., Planted/Cultivated classes from National Land 10 Cover Database 2006 (Fry et al., 2011)). Preliminary analysis indicated that random placement of multiple wells in urban and agricultural areas did not affect the streamflow depletion unless total number of wells pumping was changed. For the Haw River basin which is delineated by a total of 781 triangular elements, every single groundwater pumping well was assigned to each single triangular element as shown in Figure 1e. Thus, the amount of groundwater withdrawal for 15 public supply and irrigation were extracted from the assigned wells in the urban and agricultural areas (red and yellow colored block in Figure 1e respectively).

Fully-coupled hydrologic model -PIHM
An integrated surface water and groundwater model, PIHM, was used to estimate groundwater pumping impacts on streamflow. PIHM is a fully distributed, coupled, multiprocess model in which interception storage, overland flow, stream stage and groundwater states are solved using a semi-discrete finite volume approach (Qu and Duffy, 2007). Processes readers are referred to a supplement file (Text S1, Table S1-S2, and Figure S1). The model domain is discretized into unstructured grids (Kumar et al., 2009a). Laterally, hillslopes and rivers are discretized using triangular grids and line elements, respectively. Vertically, each groundwater layer. Soil moisture in the two unsaturated layers may vary from residual moisture to full saturation. Groundwater pumping in a given triangular grid is incorporated through a sink flux term (Kumar and Duffy, 2015). PIHM has been successfully applied at multiple scales in 15 diverse hydro-climatic regimes in both North America and Europe (e.g., Chen et al., 2015;Seo et al., 2016;Shi et al., 2013;Wang et al., 2013;Yu et al., 2013).

Hydroclimatic Data and Watershed Details
Daily precipitation and temperature data sets downloadable from Ed Maurer's webpage (http://www.engr.scu.edu/~emaurer/gridded_obs/index_gridded_obs.html) were used in this 20 study. Readers may refer to Maurer et al. (2002) for details regarding these data. These datasets are gridded observations at 1/8 degree spatial resolution and the relevant grids (total 58) Haw River basin were obtained through the USGS national map viewer and download platform webpage (http://nationalmap.gov/viewer.html/). Land cover/soil classification data were also downloadable from the geospatial data gateway of Natural Resources Conservation Service (NRCS) webpage (http://datagateway.nrcs.usda.gov/). All spatially distributed data sets were then converted into raster format for loading into PIHMgis (Bhatt et al., 2014), an open-source 20 GIS-hydrologic model framework that automatically extracts hydrogeological data sets and meteorological forcings, and maps them onto model grids.

PIHM calibration/validation
Parameters such as horizontal/vertical hydraulic conductivities, van Genuchten drainage, and porosity were calibrated with observed groundwater levels and streamflow for the Haw River basin using daily observed precipitation and temperature time series at 1/8 th degree for the calibration period from 1956 to 1980. Because of the computationally demanding nature of the 5 model, instead of using an automated calibration strategy that often requires search of the entire parameter space, manual adjustment was performed on hydro-geologic parameters such as hydraulic conductivity, macro-porosity, and soil retention parameters so that the values stay within physically reasonable range. The range of calibrated parameters are presented in Table 1.
Initial conditions on river segments (e.g., depth of the river) were given by observed data sets. 10 However, initial values of sub-surface conditions were specified due to lack of data. Given the model was calibrated for 25 years after 5 years of spin-up period (1951)(1952)(1953)(1954)(1955) and validated for the next 25 years, the specified initial conditions had a limited role in model performance ( Table   2).
The simulated monthly streamflow and groundwater level were compared with the     Apart from quantifying the differences in hydrologic variables due to pumping, we also 5 focused on how streamflow and groundwater storage recover back once the pumping stops, which we call the "no-pumping state". Since potential groundwater pumping scenarios in this study are assumed as "transient" pumping series (described in section 3.3), "no-pumping state" mentioned above is different than the "equilibrium state" discussed in previous studies (e.g., Barlow and Leake, 2012;Garner et al., 2013;Kendy and Bredehoeft, 2006;Konikow and Leake, 10 2014; Leake and Pool, 2010;Scibek et al., 2007). The term "equilibrium state", as used by  Barlow and Leake (2012) and others, meets the following criteria: water levels no longer decline in response to pumping, the cone of depression does not expand any further, and the aquifer is in a new state of equilibrium in which the pumping rate is equal to the amount of streamflow depletion. The term "no-pumping state" in this study implies that streamflow and groundwater levels are restored back to the same levels as those simulated without groundwater pumping. 5 Plus, "no-pumping state" and initial conditions (ICs) are different. ICs are initial conditions on river stage, soil moisture depth, groundwater level, etc. We define the time required for "nopumping state" to be attained after a transient pumping as the "restoration time". For example, if restoration time of streamflow (groundwater) is estimated to be 18 months, it implies streamflow (groundwater) have returned to base scenario ("no-pumping") after 18 months. We argue that

Uncertainty envelope estimation by perturbing initial conditions
Spatial variability of initial conditions can induce uncertainties on hydrologic variables especially for short-term simulation (e.g., a year simulation or less) (Moradkhani et al., 2005) Initial conditions of land surface models also influence streamflow forecasts development (Li and Sankarasubramanian., 2012;Sinha and Sankarasubramanian, 2013;Yossef et al., 2013;Li et 20 al., 2016). Given that a small difference in initial conditions such as soil moisture and groundwater storage could substantially alter the estimated groundwater time, we propose here an uncertainty estimation framework that perturbs the initial conditions under "no-pumping" 3) Uncertainty envelope on streamflow and groundwater level were then estimated as the difference between the base simulation series and each perturbed series obtained with 10 different initial conditions. Note that groundwater level is the spatial mean groundwater level over the entire watershed. We calculated spatial mean value over all elements.
Difference for each series is computed as follows. Our initial analyses on comparing the time series of hydrologic variables under 15 "pumping" and "no-pumping" showed large difference in the estimates of restoration time even though both scenarios were provided with the same climate forcings (shown in Figure 5). From times estimated using uncertainty envelope and without uncertainty envelope increased as the pumping volume increased. Given that the initial conditions of PIHM get changed once we start pumping, it is unreasonable to expect time series of various hydrologic variables under "pumping" and "no-pumping" -even though both analyses are forced with the same climate forcings -to have the same values after the system being fully recovered from a drought. The pumping" has to be tested.
We also evaluated a different noise term in (1) to make it spatially correlated by considering it as multivariate Gaussian distribution. However, it did not result in any changes in   . Demonstration of uncertainty envelope (differences between the simulation driven by "default" initial conditions and the 30 simulations driven by "perturbed" initial conditions) in estimating streamflow and groundwater depletion by comparing the "pumped" time series with the "no-pumped" time series, PRCP and GW are abbreviation of precipitation and groundwater, respectively. . Differences in estimation of restoration time for groundwater level between "considering uncertainty envelope" driven by perturbed initial conditions (blue circle) and "noconsidering uncertainty envelope" (red circle). Details on the pumping rates are described in    Figure 7 shows the annual precipitation for the three climate scenarios: the base scenario (observed forcing), the prolonged wet scenario, and the prolonged dry scenario. The logic for 15 considering these climate scenarios stems from the fact that the lag-1 autocorrelation on annual precipitation being close to zero, each year has equal probability of occurring. The autocorrelation is presented in Table S3 in a supplement file. Hence, we consider a low probable, but two realistic events of consecutive three wet years and three dry years to understand how pumping would impact streamflow depletion and groundwater storage under prolonged wet and

Groundwater pumping impacts on streamflow and groundwater depletion
BNA pumping scenarios (1 -5 BNAs -see section 3.3 for details) were imposed on the base simulation to evaluate groundwater pumping impacts on the hydrologic system. Changes in simulated variables -overland flow, infiltration, baseflow, evapotranspiration, and groundwater storage -are calculated as the difference between "pumping" and "no-pumping" conditions. We  Table 3. However, overall changes in sub-

Groundwater pumping impacts under potential climate scenario
Along with the baseline climate scenario (observed forcing), two potential extreme climate scenarios -wet scenario and dry scenario (shown in Figure 7) -were analyzed.

Restoration time
Restoration time of streamflow at outlet and groundwater level under the observed climate, wet climate, and dry climate scenarios corresponding to the different pumping volumes are presented in Figure 11a and 11b, respectively. The restoration time is defined as the time required to recover to normal conditions (i.e., within the uncertainty envelope of the baseline 5 simulation without pumping as seen in Figure 5) after the cessation of pumping.
It has been demonstrated that both streamflow and groundwater level recovered quickly (slowly) to normal conditions under wet (dry) climate in comparison to the observed climate across all pumping volumes (1 BNA to 5 BNA). The term "quickly" and "slowly" were used for relative comparison of restoration time under different climate conditions. The restoration time 10 increased with pumping volume for all three climate scenarios. The rate of increase was more pronounced for the dry scenario, moderate for the original scenario, and nearly negligible for the wet scenario ( Figure 11). Thus, prolonged dry conditions weakened the resilience of streamflow and groundwater systems due to limited water availability. Further, the relationship between restoration time and pumping volumes is non-linear while the restoration time increased 15 monotonically for large pumping volumes. It is important to mention that over-exploitation of groundwater pumping can also impact hydro-ecosystem since it would take long time for the inflow and groundwater to come back to normal condition. This threat of over-exploitation to groundwater resources actually has been a serious issue in the regions affected by dry climate regime (e.g., Southern California and Arizona) (Draper, 2015;Lustgarten, 2015). Thus, any 20 efforts in pumping during drier conditions should also consider potential streamflow depletion and longer restoration time of groundwater level. We intend to consider this as part of our future research.

Discussion and Concluding Remarks
Since transient pumping impacts at the river basin scale could be estimated only through modeling efforts, we argue that it is important to consider relevant uncertainties. Hence, we proposed the uncertainty framework that rigorously incorporates uncertainties in initial 5 conditions and input variables (Li et al., 2009;Sinha and Sankarasubramanian, 2013;Yossef et al., 2013). We proposed an uncertainty framework that estimates the restoration time of groundwater and surface water systems by comparing the distribution of streamflow under "nopumping" with the streamflow obtained under "pumping". We demonstrated the application of this framework in estimating the restoration times and the associated depletion volumes due to 10 groundwater pumping for a humid basin in NC by considering different climatic and demand conditions. Using the uncertainty framework, the resilience of the system was examined along with the changes in water-budget components due to pumping using the integrated surface-water and groundwater model -PIHM. Pumping impacts under potential climate scenarios were analyzed to understand how different sequences of wet and dry years could affect the restoration 15 time and depletion volumes of streamflow and groundwater levels. We first assessed how small differences in initial watershed conditions result in significant overestimation of restoration times. Hence, instead of simply comparing two time-series of streamflow and groundwater developed under "pumping" and "no-pumping", we proposed an uncertainty framework that perturbs initial conditions to obtain a null distribution of streamflow and groundwater under "no- 20 pumping" for comparing the time series of water-budget components obtained under "pumping".
We also suggest consideration of additional sources of uncertainty such as model, parameter and input uncertainties for developing robust estimates of restoration time (Li and Sankarasubramanian, 2012;Li et al., 2016).
Our analyses show that the restoration time of groundwater is more sensitive to climatic conditions (wet vs dry) as opposed to the pumping rates. The restoration time of groundwater and streamflow changed little with continued wet conditions. As expected, the restoration time 5 significantly increased with consecutive dry years. Analyses of water-budget components and state variables revealed several critical insights. Evapotranspiration from the watershed reduced due to low soil moisture (compared to non-pumping conditions) under pumping. Similarly, baseflow contribution also reduced due to reduced groundwater storage under pumping over the watershed. Since the water table is shallow (the spatial mean depth to water table was 10 approximately 2 meters in the model), groundwater pumping influences reductions in evapotranspiration, baseflow and overland flow. After the cessation of pumping, increased replenishment of groundwater storage due to enhanced infiltration also reduced overland flows.
Thus, an integrated assessment of watershed conditions due to groundwater pumping was performed to understand how the hydrologic variables and storages of the PIHM changed due to 15 pumping.
This study emphasized the importance of analyzing potential pumping impacts on watershed particularly on the ability of surface water and groundwater systems to recover during a drought. We found significant uncertainty in recovery time exists if uncertainty in initial conditions is not properly considered. Lack of reliable groundwater pumping data sets and  (Devineni and Sankarasubramanian, 2010a, b). Thus, for future study, any 5 efforts to manage surface water and groundwater resources during drought conditions should consider conjunctive management by quantifying the watershed resiliency conditioned on climate forecasts (Devineni and Sankarasubramanian, 2010a, b).
Analyses for the humid basin in a rainfall-runoff regime also provide insights on potential implications for other regimes. In arid settings, it is natural to expect the competition for 10 evapotranspiration to increase due to increased potential energy availability (Sankarasubramanian and Vogel, 2005). Hence, groundwater pumping in arid settings is further expected to increase the restoration time as the potential for evapotranspiration increases. This in turn is expected to increase the streamflow and groundwater depletion. However, in snowmelt dominated regimes in humid basin, the opportunity for infiltration increases. Hence, we can 15 expect relatively a quicker recovery time after pumping. These are interesting aspects that require further research to draw a synthesis on how groundwater pumping impacts various waterbudget components and watershed resiliency over the Coterminous US (CONUS). We intend to consider this as part of our ongoing research in understanding the human and climate influence watershed resiliency over the CONUS.