Water Balance to Recharge Calculation: Implications for Watershed Management Using Systems Dynamics Approach

Groundwater depletion in the face of growth is a well-known problem, particularly in those areas that have grown to become dependent on a declining resource. This research comprises a broad synthesis of existing water resources data, to understand the long-term implications of continued growth in water demand on groundwater dominant water resources, and to develop a tool for sustainable water management. The Palouse region of Washington and Idaho, USA. (approximately 60,000 people in a rural setting) is entirely dependent on groundwater from two basalt aquifers for potable water. Using the systems dynamics approach and a water balance that considered the entire hydrologic cycle, a hydrologic model of these aquifers was developed, tested and applied to simulate their behavior over a 150 year time period assuming the current infrastructure does not change. With 1% population growth and current water extraction rates, the results indicated the upper aquifer use may be sustainable, while the lower aquifer use is likely unsustainable in the long term. This study also shows that uncertainties in key aspects of the system create limitations to groundwater management.


Introduction
The continuing extraction of groundwater from the aquifers for drinking, agricultural and industrial usage is contributing to groundwater depletion in many parts of the United States and the world [1][2][3][4][5][6][7][8][9][10][11].Managing groundwater with limited information is complicated, and dealing with uncertainties not only requires significant scientific effort, but can lead to ineffective management that is slow to respond to observed depletion.Information is generally limited because of complex hydrogeology, difficulties in and costs of data collection, and the need for frequent monitoring.For example, recharge is considered to be one of the most difficult components of a water-balance to quantify accurately as it varies significantly in both temporal and spatial dimensions [12][13][14][15][16][17][18][19].Effective management of groundwater aquifers is heavily influenced by the accurate estimation of recharge.Another, perhaps under-appreciated characteristic of large confined aquifers that affects long-term management is storativity, as it defines the volume of water that is or could be stored for future use.Together, recharge and storativity and their uncertainty are critical for groundwater management.
The Palouse region is a semi-arid area located along the border of northern Idaho and eastern Washington, is solely dependent on groundwater for drinking water.Two distinct confined groundwater aquifer systems exist in the Palouse Basin (PB), the upper Wanapum (WP), and the lower Grande Ronde (GR).The GR aquifer supplies the potable water for approximately 70 percent of Moscow, Idaho and 100 percent of Pullman, Washington.The WP aquifer is on average 110 m below ground surface and is thought to receive some recharge, while the GR aquifer is on average approximately 290 m below ground surface with indeterminate recharge rates.Since the first well was installed over 100 years ago, water levels in the GR have been decreasing 0.3-0.45m annually indicating long-term unsustainable use [20].Much effort has been expended to understand basin hydrogeology and groundwater recharge of the PB (e.g., [21][22][23][24][25][26][27][28][29][30][31]).To estimate basin-scale recharge of the PB, Dhungel [29], conducted a water balance at the land surface based on 30-year average precipitation, streamflow, and evapotranspiration; it is assumed that this recharge is received by the upper WP.The lower GR is often thought to receive little to no recharge.This assumption is supported by some studies ( [28] and [32] but many believe the existence of recharge to GR ( [26,[33][34][35]).
In the 1950s, the systems dynamics method was formally created by J. W. Forrester at Massachusetts Institute of Technology [36], which has proved to be useful to simulating any system that can be thought of in terms of stocks and flows.The approach is often used for shared vision planning of natural resources, a process which involves discussion and debates as part of the comprehensive decision making process [37].The systems dynamics approach facilitates understanding of the behavior of complex systems over time.As such, the use of systems dynamics in water resources planning accelerated in the 1990s, for example in water resource planning [38][39][40][41][42][43][44][45][46][47][48][49] and for groundwater management [50][51][52].Mirchi et al. [53] synthesized the integrated water resources modeling efforts.Dhungel [29] and Beall et al. [45] discussed water resources management efforts in the Palouse region using a systems approach.
Thus, the overall objective of this study was to synthesize the currently available hydrologic data, and to incorporate the results of this synthesis as a lumped water balance model using a systems dynamics approach.This model is then used to simulate future water resources for various growth scenarios and to show the effects of uncertainties.

Study Area and Data
The majority of land in the Palouse Basin is used for dryland agriculture (grains and legumes) and forest.The soils are primarily silty loam loess, wind-transported soils that form rolling hills, which overlies Columbia River basalt flows.Palouse loess elevations range from 1010 m along the eastmost boundary to 520 m in the west.The climate is semi-arid, with precipitation falling primarily in the winter months, and the summers are relatively hot and dry.The total population of Palouse region was about 51,000 (2006) which is comprised of small cities like Moscow, Potlatch, Pullman, Viola, Potlatch, and Colfax.Approximately 9.5 percent of the total water demand of the PB was fulfilled by the WP and rest by the GR.Total water extraction from the PB was about 2.77 billion gallons in the year 2005 [54].Water extraction from the GR (2.51 billion gallons/year) was approximately nine times greater than the WP (0.26 billion gallons per year) in the year 2005.The per capita per day water use of the PB is approximately about 150-180 gallons [20,21], but recent studies show a decline in water use, which is about 120 gallons [55] that is used in the simulation.
Population data are acquired from the United States Census Bureau and city sources.The population of this area is estimated to grow by 1% annually [20].A longer simulation period, i.e., 150 years were chosen in contrary of 50-100 years because it is assumed that the population growth of these rural settings would rise constantly rather than a large spike.Figure 1 shows the location map of the study area and groundwater and surface area regions.The groundwater aquifers of PB are located within the basaltic Columbia River flows.Water level of principal well in the basin is about 686 m (above mean sea level) in 1990s [56].Even though, there is water below 488 m (i.e., 1600 ft) but uncertainty of water volume increases significantly [45].The existence of the upper WP aquifer seems significant in all groundwater basins with comparatively thin layer in Pullman (46 m) [57].The Moscow WP is productive for groundwater extraction whereas the Pullman WP is unproductive [58].The uppermost layer of the Palouse Basin is composed of loess which is basically a deposit of wind-blown silt.The composition of these aquifers is more than 60 percent basalt, with the rest being sediments including silt, clay and sand.According to the geographic variations, the groundwater regions are divided into Palouse, Colfax, Viola, Pullman, Moscow and Uniontown regions (Figure 1, [57]).The radiocarbon age of GR groundwater basin varies between 8500 and 31,000 years apparent Radiocarbon Age (BP) [32], between 13,000 and 26,400 years (BP) for the lower portion of the GR and 4400 to 11,800 (BP) for upper portion of the GR [59] and approximately 20,000 years [45].regions are divided into Palouse, Colfax, Viola, Pullman, Moscow and Uniontown regions (Figure 1, [57]).The radiocarbon age of GR groundwater basin varies between 8500 and 31,000 years apparent Radiocarbon Age (BP) [32], between 13,000 and 26,400 years (BP) for the lower portion of the GR and 4400 to 11,800 (BP) for upper portion of the GR [59] and approximately 20,000 years [45].The question of the storativity of the basaltic rock of PB aquifer and Columbia basin is a particular interest of geologist, water managers and hydrologist.The two closely located state universities of Idaho and Washington along with other entities are working to better understand the aquifer properties and hydrogeology of the basin.Studies show that storativity of PB aquifer basalt ranges approximately between 10 ´3 and 10 ´5 based on aquifer discharge tests [60], 1.5 ˆ10 ´4 ˘0.2 ˆ10 ´4 [61], 3 ˆ10 ´5 to 3 ˆ10 ´4 of GR [62] and 0.67 ˆ10 ´3 of WP [63].Moran [62] had compiled the storativity estimates from previous PB investigations.Even though the storativity values of the aquifer are vital because it governs the initial volume and overall life of aquifers, recharge to the aquifers plays a critical role.The short term drawdown and recovery  of WP indicates that recharge was the prevailing factor of the dynamics of WP.The aquifer recharging process can be less evident in GR because of the larger surface area, complex hydrology, and hydraulic and hydrologic connectivity among the other groundwater regions, but the inconsistency between the drawdown and extraction is indicating that GR might be actively recharging.Analyzing initial volume of aquifers, the current decline rate the aquifer and the long term pumpage, we utilized the storativity value of 10 ´3-10 ´4 [61] for the modeling purposes.Smaller storativity values imply less water stored in the specific saturated depth of aquifer as well as only a little amount of water that can be extracted from drawdown.
To conduct basin scale water balance, different sets of climatic and surface water data is needed.For calculating evapotranspiration, land use maps, soil maps, elevation maps and PRISM maps [64] were used.To compute the surface water, daily discharge from USGS gauging stations was used from 1970 to 2000.A watershed was delineated using 10 m resolution Digital Elevation Map (DEM) of PB.Colfax (i.e., USGS gauging station 13346100, Palouse River at Colfax, WA) was taken as the downstream point for delineating watershed.Five more USGS gauging stations were used for delineating sub-watersheds with a total area of 2044 km 2 representing entire PB (Table 1).The basin was divided into six surface water sub-basins defined by the United State Geological Survey (USGS) surface water gauging station locations (Figure 1).

Methodology
In this section, discussions are carried out of the components of the water balance, governing equations and the applied boundary conditions.

Water Balance Analysis
The hydrology section is divided into surface water and groundwater.Surface water hydrology is described by mean annual precipitation, surface runoff and evapotranspiration; recharge is computed by water-balance, assuming that the average soil moisture storage change in the unsaturated zone is zero.Equation (1) shows a general form of the water-balance approach used to calculate recharge to WP [65][66][67][68].The recharge was calculated using the basic water-balance approach over the 1971 to 2000 time period (2044 km 2 ).Scanlon et al. [12] classified recharges broadly as physical based and tracer or numerical modeling approach.Scanlon et al. [69] categorized the techniques for estimating recharge for surface water, unsaturated zones, and saturated zones.Risser et al. [70] used unsaturated-zone drainage collected in gravity lysimeters, daily water balance, water-table fluctuations in wells, and equations of Rorabaugh for estimating ground-water recharge at a humid-continental climate in east-central Pennsylvania.Lee et al. [71] discussed the standard techniques for regional groundwater recharge calculation based on water-balance model and parameter-value adjustment of groundwater flow models.Simmers [72] defined recharge as short-term, seasonal, perennial and historical based on the different time scales.A water budget or balance approach is the largest class techniques to estimate recharge [73].
The term P (L 3 /T) represents precipitation, ET (L 3 /T) is the evapotranspiration, Q s (L 3 /T) is the discharge of the river (surface water inflow), ∆S (L 3 /T) is storage change in the unsaturated zone, and R (L 3 /T) is the recharge.
Mean areal annual precipitation in Equation ( 1) was computed from Parameter-elevation Regressions on Independent Slopes Model (PRISM) maps of girded data with 30-arcsec (800 m) NORMALS .The estimation of evapotranspiration in equation 1 followed the prediction methodology by the Thornthwaite and J.R. Mather approach.Thornthwaite-Mather [74] is a lumped model where the entire watershed is treated as a single unit and soil water status is tracked through time.Specific parameters used in this method are rooting depth, available soil water storage depth, crop coefficient and maximum canopy storage amount.The potential evapotranspiration was calculated using the Hargreaves approach [75].Daily discharge data  from the USGS gauging stations were used to calculate surface runoff (Equation ( 1)).Surface water watersheds are delineated solely from USGS gauging station locations in order to perform mass-balance computations.Table 1 also includes the USGS gauging stations, period of data availability and the linear regression equations that are used to estimate missing data such that the water balance is based on a consistent period of record.
The groundwater volume of the aquifers (V (L 3 )) were estimated from potential groundwater drawdown of average saturated depth (H (L)), surface area of the aquifer (A (L 2 )) and storativity (S) [76].
Water-balance of WP (Equation ( 3)) is computed by balancing the initial volume of WP (V WP , Equation ( 2)), recharge to WP (R WP ), discharge from WP (D WP ) that does not reach to GR, WP pumpage (P WP ), and recharge to GR (R GR ).It is assumed that about 12% of water is extracted from Wanapum and rest from GR.
Change in storage in WP (∆S WP ) can be written as of Equation ( 4).
V WP ptq ´VWP pt ´dtq " ∆S WP " 0 pA f ter 1990q In Equation ( 3), all the components are known except D WP and R GR .Discharge from Wanapum (D WP ) that does not reach to GR is computed from the following equation simplifying Equation ( 3) and (4) (Equation ( 5)).
Water-balance of GR is computed by balancing the initial volume of GR (V GR ), recharge to GR (R GR ), and GR pumpage (P GR ).In Equation ( 6), all the components are known except R GR .
Change in storage in GR (∆S GR ) can be written as of Equation (7).The term ∆h GR [L] represents drawdown of GR and A GR is surface area of GR.
Recharge to GR can be computed simplifying Equations ( 6) and ( 7) (Equation ( 8)).With this assumption, any discrepancy in estimation of storativity or surface area of GR can affect the overall accuracy of the estimated recharge.Similar approach (Equation ( 8)) had been utilized to test how the system parameters affect the recharge estimates in GR [61,62].The maximum recharge to GR (R GR (max) ) is limited to 4 cm/year not to surpass the recharge from water balance.
Once the recharge reaches to the maximum value, drawdown can be back calculated using following equation.
The simulation process is stopped when GR aquifer no more able to fulfill the entire demand.

Modeling Approach
A schematic of the PB watershed is shown in Figure 2a where the surface area of the PB, the productive WP region and GR region varies.These surface area of aquifers is used to compute the initial volume of groundwater, which represents the start of the simulation period which is 1930.Prior to 1930, it is assumed that there was no significant overdraft of groundwater in both aquifers.As it is generally thought that GR get recharged from the fractures, so the actual pathways of recharge to GR is insignificant.GR regions are assumed as non-leaking reservoirs.Studies show the GR aquifer system of the Moscow, Pullman, Colfax and Garfield are hydraulically connected [77].Based on this finding, we combined the groundwater regions of GR to develop a lumped GR hydrological unit (728 km 2 , (610 km 2 , suggested by [59]) for water storage and water balance purpose.Similarly, lumped WP region is associated with Moscow WP (82 km 2 ) assuming that the rest of WP groundwater regions are unproductive (Figure 2a).If the surface area of WP and GR are assumed to be identical of PB (2044 km 2 ), the groundwater storage and extraction (Figure 2b) will be increased.Simulations will be carried out how these affect the final results.
Figure 3a shows the causal loop diagram of the Palouse basin watershed in its simplest form.An increase in precipitation will increase recharge to WP and GR.Increased population in PB would put pressure to both WP and GR which will ultimately reduce the storage of aquifers by further increasing the drawdown.
A systems model of the PB watershed using STELLA software [78] shown in Figure 3b.Historical and current observations show (1930 to current) that the average annual drawdown of GR is about 0.37 m irrespective of approximately five times population growth (10,000 to 50,000).The annual water level drop in GR actually slightly decreased from 0.45 m to 0.3 m in 1990s.These water level decline trends are indicating that multiple factors (probably recharge, connectivity of groundwater regions, storativity, porosity, resident time, groundwater dating etc.) are influencing along with the population growth.Beall et al. [45] purposed dual porosity to estimate the aquifer decline rate of GR using resident time concept.Dhungel [29] had conducted an extensive sensitivity analysis to understand the aquifer behavior of PB varying storativity, surface area and recharge, but because of the uncertainty in recharge and storativity values, the effort provided a little insight to understand the life of these aquifers.In this study, we utilized the indirect substitution method rather direct sensitivity analysis to estimate recharge to GR (Equation ( 8)).According to Equation ( 8), GR pumpage constantly exploits WP in the form of recharge to balance the deficit between the GR pumpage and drawdown.
productive WP region and GR region varies.These surface area of aquifers is used to compute the initial volume of groundwater, which represents the start of the simulation period which is 1930.Prior to 1930, it is assumed that there was no significant overdraft of groundwater in both aquifers.As it is generally thought that GR get recharged from the fractures, so the actual pathways of recharge to GR is insignificant.GR regions are assumed as non-leaking reservoirs.Studies show the GR aquifer system of the Moscow, Pullman, Colfax and Garfield are hydraulically connected [77].Based on this finding, we combined the groundwater regions of GR to develop a lumped GR hydrological unit (728 km 2 , (610 km 2 , suggested by [59]) for water storage and water balance purpose.Similarly, lumped WP region is associated with Moscow WP (82 km 2 ) assuming that the rest of WP groundwater regions are unproductive (Figure 2a).If the surface area of WP and GR are assumed to be identical of PB (2044 km 2 ), the groundwater storage and extraction (Figure 2b) will be increased.Simulations will be carried out how these affect the final results.Figure 3a shows the causal loop diagram of the Palouse basin watershed in its simplest form.An increase in precipitation will increase recharge to WP and GR.Increased population in PB would put pressure to both WP and GR which will ultimately reduce the storage of aquifers by further increasing the drawdown.We assumed that the current water depletion trend will continue for a certain period of time in the future, eliminating one of the unknowns in Equation ( 6), which ultimately helps to calculate recharge to GR.The vital question is if this trend continues, then how long it last and how will it respond to the increased population growth.Obviously, increased population growth will accelerate the groundwater extraction and shorten the life of these aquifers, however, the water decline trend of GR may not be linearly proportional to population growth as we believe there might be active recharge.Our goal here is to develop a model which can sustain the current water decline rate for the existing population growth (1%) as well as the increased growth.The only way it can be achieved is from additional contributions from WP to GR as recharge, but recharge to GR from WP is the function of water balance of PB.The simulation will find out the cutoff point where these water balance equations and the hypothesis are not further applicable.The actual water level trend in the simulation may be less intuitive, but these cutoff points will indicate the optimum duration that this trend can exist with these hypothesis and assumptions.The model will face two types of threshold or cutoff points a) basically when the demand of GR i.e., GR pumpage is higher than GR storage (deficit threshold) b) when the recharge to GR reach the maximum allowed values (recharge threshold).A systems model of the PB watershed using STELLA software [78] shown in Figure 3b.Historical and current observations show (1930 to current) that the average annual drawdown of GR is about 0.37 m irrespective of approximately five times population growth (10,000 to 50,000).The annual water level drop in GR actually slightly decreased from 0.45 m to 0.3 m in 1990s.These water level decline trends are indicating that multiple factors (probably recharge, connectivity of groundwater regions, storativity, porosity, resident time, groundwater dating etc.) are influencing along with the population growth.Beall et al. [45] purposed dual porosity to estimate the aquifer decline rate of GR using resident time concept.Dhungel [29] had conducted an extensive sensitivity analysis to understand the aquifer behavior of PB varying storativity, surface area and recharge, but because of the uncertainty in recharge and storativity values, the effort provided a little insight to understand the life of these aquifers.In this study, we utilized the indirect substitution method rather direct sensitivity analysis to estimate recharge to GR (Equation ( 8)).According to Equation ( 8), GR  The drawdown of the groundwater aquifers can vary both spatially and temporally and needs an enormous number of monitoring wells to understand basin scale drawdown.As a part of model validation, we mimicked drawdown behavior of these aquifers based on the historical and current observations.To replicate the drawdown behavior of the WP between 1930 and 1960, a constant drawdown (∆h WP : approximately 0.1 m/year for 30 years) was adopted.During this period, the majority of water was extracted from the upper WP which generated rapid drawdown in WP.When the extraction from WP was dwindled to overcome this drawdown, the water level of WP started attaining back to an earlier state at 1990s [45].During this period (1960-1990), increasing ∆h WP was adopted to balance the previous loss.Finally, from 1990s, ∆h WP is considered to be zero as WP is nearly constant.To replicate the GR decline trend, a constant value of ∆h GR (i.e., 0.45 m/year) was adopted between 1930 and 1990 and after 1990, 0.3 m/year.

Water-Balance and Recharge to Aquifers
In this section, discussion is carried out of individual component of the water balance of particular sub-watersheds.As we are generating a lumped model, there may little importance of particular components of sub-watersheds, but it is intended to demonstrate how components of water balance can vary in a fairly small watershed.The maximum precipitation value within the watershed is approximately 85 cm in the Palouse River sub-basin above Potlatch, Idaho, while the lowest value of approximately 59 cm is observed at South Fork above Colfax, Washington.The mean areal precipitation over the entire watershed is approximately 71 cm.The maximum mean areal surface runoff is approximately 29 cm in the Palouse River near Potlatch, Idaho, while the lowest value of 5 cm is observed in the Palouse River at Palouse, Washington.The mean areal runoff over the entire watershed is approximately 17 cm.The maximum mean areal ET is approximately 54 cm in the Palouse River near Potlatch, ID, while the lowest value of approximately 45 cm is observed in the South Fork above Colfax, Washington (North Fork).The mean areal ET over the entire watershed is approximately 49 cm.Table 2 shows the estimated surface runoff, evapotranspiration and precipitation of each sub-watershed of WP aquifer.The areal mean precipitation, evapotranspiration and runoff of entire watershed was 71, 49 and 17 cm resulting recharge of 4.7 cm.The wide variation in calculated recharge from the surface watersheds (Table 2) suggests that replenishment of the underlying WP aquifer would also be spatially variable.Some portions of the WP aquifer will be more easily recharged than others due to differences in recharge rate.To understand these variations, Dhungel [29] developed a hydrologically separated model based the geological separation of groundwater regions.A water-balance at the land surface was used to estimate recharge to the Wanapum by spatially sub-dividing surface water areas and groundwater regions.However, the population center of PB lies within 7 miles i.e., Moscow and Pullman area and the hydrogeological investigations indicated that Moscow WP is only productive, the hydrologically separated model faces similar constraints to the lumped models [29].Because of significant heterogeneity in the GR, traditional numerical model of groundwater flow (e.g., MODFLOW) had offered little insights into the overall nature of the aquifer [22,45].
Table 3 shows some of the documented WP recharge values based on various approaches since 1960s.These studies suggest the recharge rates of upper WP varies from 0.3 to 19.7 cm/year.The large differences in these recharge calculations increase uncertainty in water storage in WP.Based on the expert elicitation and Bayesian Model, Reeves [30] found the estimated WP areal average recharge to be 5.1 cm¨year ´1 with an uncertainty of ˘4.6 cm¨year ´1.Based on the recharge calculation (Table 2), the water-balance shows 383 billion gallons (BG) of water are contributed by precipitation, 93 billion gallons of water are converted into surface runoff, and evapotranspiration accounts for 267 billion gallons, resulting in 22 billion gallons of recharge to the WP for the year 2005.This amount of recharge (approximately 4 centimeters per year) corresponds to an average rate over the entire area of 2044 kilometers.This estimated recharge amount to WP aquifer is approximately eight times larger than the present water extraction of entire PB (2.77 billion gallons, 2005).This may be the possible indication, in part (1) the computed recharge is likely overestimated (2) entire recharged water does not reach to WP (3) GR get constant recharge from WP (4) possible inconsistency of surface areas of PB, WP and GR.The computed recharge from the water-balance is in the range of the previously documented values (Table 3).As per the assumptions, the portion of recharged water might be discharged from WP while the rest might be possibly recharging to GR (Equations ( 3)-( 5)).

Aquifer Storage
One of the frequently raised questions while managing the groundwater is the aquifer storage capacity and possible life of the aquifers.However, the larger uncertainty in storativity and hydrogeology can put these topics in shadow and the usefulness of these simple models can also be questioned.The estimated total initial volume of the WP and GR is about 0.3 billion gallons (surface area about 82 km 2 ) and about 5 billion gallons (surface area of 728 km 2 ), respectively, when potential drawdown at the bottom of the aquifers with storativity of 10 ´4 (Table 4).This total volume of groundwater in the WP and GR represents the maximum water in the aquifers irrespective to the present water level.The average maximum potential drawdown of the Moscow WP is about 137 m, Moscow GR is 271 m and Pullman GR is 305 m at the bottom of aquifers.Now, coming back to the original question if these estimated aquifer storage capacity is reasonable, while comparing the actual water extraction.Water extraction from these aquifers far outreaches the estimated storage capacity or our understanding to storativity and recharge needs extensive revisiting.As an optimistic scenario, simulation was also carried out of the storativity of 10 ´3 (based on the published storativity values).The higher storativity values not only imply a larger storage capacity, but also a larger amount of water literally extracting from the current drawdown.Similarly, smaller storativity indicates a large amount groundwater needs to be replenished as recharge to balance the current pumpage.The amount of water that is continuously declining from GR is about 0.06 billion gallons (storativity = 0.001, 0.3 m/year) and 0.01 billion gallons (storativity = 0.0001, 0.3 m/year), which is a tiny fraction of the current extraction.It indicates the substantial contribution of recharge in both scenarios.

Simulation of Groundwater Volume of the Aquifers
Based on the water balance (Equations ( 1) to ( 9)), simulations were carried out to find the response of these aquifers.Figure 4a,b show the simulation results of groundwater storage change between 1930 and 2150 with storativity values of 0.0001 and 0.001, respectively.The population of the PB was validated in 2005 and 1%-2% increase in population was adopted after that.The percentage value in the parenthesis in the figure shows the applied population growth in the simulation.As expected, simulation result shows that the WP decreased until 1960s and reverts back its original level in the 1990s.After 1990s, the water level of WP is constant as the model considers no drawdown.The aquifer dynamics of upper WP may have little interest compared to the lower GR as it is fairly constant, but the entire recharge process to GR is dominated by WP.In the model, GR adopts a constant decline of 0.45 m/year until 1990s and after that, a decline about 0.3 m/year until the maximum allowed recharge reached, i.e., 4 cm/year.Figure 4a (smaller storativity value of 0.0001) shows the validity of these assumptions are until about 2070 and 2040 to 1% and 2% population growth, respectively.Because of the smaller storage capacity of GR, i.e., 5 billion gallons, it will not be able to fulfill the entire demand after these cut off points.These cutoff points do not imply that the groundwater cannot be further extracted, but probably only a portion of demand can fulfilled or alternative sources need be sought out to fulfill the entire demand.When a 2.5% population growth is applied, there was a significant decline in water level of GR after 2050 based on the analysis conducted by Beall [45].
In an optimistic scenario (with larger storativity i.e., 0.001), these assumptions can be valid until 2090 for 2% population growth and beyond 2150 to 1% population growth.The actual shape of these simulations might be less intuitive as we have adopted constant drawdown in the simulation, but looking back 100 years drawdown, this assumption can be acceptable for near future.Figure 4b shows the illustration of the curve fitting of aquifer decline if constant drawdown replaced by variable drawdown ( for 1%, 2% population growth or somewhat between all these simulations, purple line) assuming the cutoff points is reasonable.With a smaller storativity value, to maintain the current pumpage from GR, either larger groundwater regions or larger recharge rate is needed [62].After analyzing the results (Figure 4a,b), it is not pragmatic to increase the population growth lager than 2%.Sensitivity analysis as well as the uncertainty analysis can be conducted altering numerous parameters of the aquifer (surface area, storativity, groundwater depth, recharge, population growth etc.), but we believe that the model had sufficiently captured the behavior of the aquifer with a range of applied storativity, variable recharge rate as well as surface area.Figure 4c shows the observed long-term monthly hydrograph of GR wells at Washington State University (WSU) and GR test well (WTEST).
As expected, simulation result shows that the WP decreased until 1960s and reverts back its original level in the 1990s.After 1990s, the water level of WP is constant as the model considers no drawdown.The aquifer dynamics of upper WP may have little interest compared to the lower GR as it is fairly constant, but the entire recharge process to GR is dominated by WP.In the model, GR adopts a constant decline of 0.45 m/year until 1990s and after that, a decline about 0.3 m/year until the maximum allowed recharge reached, i.e., 4 cm/year.Figure 4a (smaller storativity value of 0.0001) shows the validity of these assumptions are until about 2070 and 2040 to 1% and 2% population growth, respectively.Because of the smaller storage capacity of GR, i.e., 5 billion gallons, it will not be able to fulfill the entire demand after these cut off points.These cutoff points do not imply that the groundwater cannot be further extracted, but probably only a portion of demand can fulfilled or alternative sources need be sought out to fulfill the entire demand.When a 2.5% population growth is applied, there was a significant decline in water level of GR after 2050 based on the analysis conducted by Beall [45].In an optimistic scenario (with larger storativity i.e., 0.001), these assumptions can be valid until 2090 for 2% population growth and beyond 2150 to 1% population growth.The actual shape of these simulations might be less intuitive as we have adopted constant drawdown in the simulation, but looking back 100 years drawdown, this assumption can be acceptable for near future.Figure 4b shows the illustration of the curve fitting of aquifer decline if constant drawdown replaced by variable drawdown ( for 1%, 2% population growth or somewhat between all these simulations, purple line) Figure 5 shows the simulation results of the possible contribution of recharge to GR from WP in centimeters.Historical and current observations show that GR is declining nearly at the constant rate even though there is a significant surge in GR extraction (from 1930 ´0.14 billion gallons to 2005 ´2.5 billion gallons).Water balance indicates that the WP may receive a considerable water from the entire PB (up to 22 billion gallons/year) compared to the tiny amount of water extraction from WP (less than a billion gallons/year), but the current water level of WP is fairly constant.Similarly, the current pumpage from GR (>2 billion gallons) largely exceeds the observed drawdown (~0.02-0.09 billion gallons) which also suggests that the GR may be actively recharging.Figure 5 shows that the majority of the recharged water from WP is either departing or potentially unusable when there is less demand in GR.However, when the demand in GR increased, more water from WP is infiltrating to GR, indicating a variable recharge rate.In a separate simulation, we tested a hypothesis of no recharge to GR which makes the GR as a sole function of GR pumpage (Figure 6).If assumed no recharge to GR, water extraction from GR should be equal to drawdown of GR.Without an acceptable storativity value, we can argue the fruitfulness of this simulation, but it will help to distinguish the role of recharge.Simulation result shows that entire groundwater from GR should have been emptied in around the 2000s (with higher storativity 0.001, Figure 6), with the current applied drawdown (0.3 m/year) and 1% population growth.Analyzing this scenario, we attempted to increase life of GR by increasing surface area of GR assuming identical to the entire PB (2044 km 2 ).This attempt does not help to increase the life of the GR because of the lack of continuous source to replenish the pumpage.This probably implies that a significant amount of water might be entering as recharge to GR to balance the deficit.In a recent study, Piersol and Sprenke [61] indicated that about 93.5% ± 2.6% of water in GR is contributed by the recharge which agrees the finding this study (Figure 7).
Figure 7 shows the contribution of recharge to GR as a percentage.It was computed based on the water that is extracted from the GR to recharged to GR. Results show that up to 96% of GR demand might had been fulfilled by the recharge between 1930 and 2005.In the future, more water is needed to fulfill the increased demand of GR to maintain the current decline rate of 0.3 m/year or GR will face a larger decline (Figures 4 and 5). Figure 7b shows that the % contribution of recharge decreased when recharge value reached the threshold value, i.e., 4 cm/year.While discussing the recharge to GR, an important concern is how these estimated recharge values comply with the age dating of GR groundwater.As mentioned earlier, some of the studies dates water from GR as early as 30,000 years (BP).The answer might not be straightforward and probably the most up-to-date explanation can be considered to be the approach adopted by [45].This study keeps on the paradigm that GR may be actively recharging based on the water balance model which is supported by some of the most current studies [61,62].Figure 5 shows that with an increase of water demand in GR, recharge from WP to GR is also consistently increasing until it reaches the threshold.Figure 5a shows that after 2060, GR no more able to fulfill the entire water demand of PB for 1% population growth and even shorter period for 2% growth, which is an example of deficit threshold implied in the model.The maximum recharge value is about about 2 cm/year within these periods.
In an instance of a higher storativity with larger storage capacity (0.001, Figure 5b), the GR will be able to fulfill the entire water demand for 1% population growth until 2140 with the assumed constant drawdown.When the recharge from WP to GR reaches the recharge threshold (Figure 5a,b), there will be a sharp decline of groundwater beyond this point (drawdown is shown by the red line in secondary vertical-axis).Beall et al. [45] utilized the dual storage and resident time concept to understand the behavior of the GR aquifer.Results confirmed that in order to maintain the current trend of drawdown in the future, either more water from WP needs to enter in GR or GR will observe larger drawdown (red line at secondary axis) if recharge to GR is constrained by other hydrogeological conditions.As discussed earlier, it is difficult to predict the actual behavior and the possible trend of groundwater depletion in the future because it might be the function of a multiple parameters, but the larger overdraft of groundwater would obviously increase the drawdown.
In a separate simulation, we tested a hypothesis of no recharge to GR which makes the GR as a sole function of GR pumpage (Figure 6).If assumed no recharge to GR, water extraction from GR should be equal to drawdown of GR.Without an acceptable storativity value, we can argue the fruitfulness of this simulation, but it will help to distinguish the role of recharge.Simulation result shows that entire groundwater from GR should have been emptied in around the 2000s (with higher storativity 0.001, Figure 6), with the current applied drawdown (0.3 m/year) and 1% population growth.Analyzing this scenario, we attempted to increase life of GR by increasing surface area of GR assuming identical to the entire PB (2044 km 2 ).This attempt does not help to increase the life of the GR because of the lack of continuous source to replenish the pumpage.This probably implies that a significant amount of water might be entering as recharge to GR to balance the deficit.In a recent study, Piersol and Sprenke [61] indicated that about 93.5% ˘2.6% of water in GR is contributed by the recharge which agrees the finding this study (Figure 7).  Figure 8 shows the comparison among the different pumping entities and recharge to GR.The GR pumpage closely followed to the recharge to GR.It is because of the dependence of pumpage to recharge and it this will increase in the future as the demand increases.At the end of the simulation period i.e., at 2150, the demand of PB (combined WP and GR) will be about 10 billion gallons (1% growth).Figure 8 shows the comparison among the different pumping entities and recharge to GR.The GR pumpage closely followed to the recharge to GR.It is because of the dependence of pumpage to recharge and it this will increase in the future as the demand increases.At the end of the simulation period i.e., at 2150, the demand of PB (combined WP and GR) will be about 10 billion gallons (1% growth).Figure 7 shows the contribution of recharge to GR as a percentage.It was computed based on the water that is extracted from the GR to recharged to GR. Results show that up to 96% of GR demand might had been fulfilled by the recharge between 1930 and 2005.In the future, more water is needed to fulfill the increased demand of GR to maintain the current decline rate of 0.3 m/year or GR will face a larger decline (Figures 4 and 5). Figure 7b shows that the % contribution of recharge decreased when recharge value reached the threshold value, i.e., 4 cm/year.
Figure 8 shows the comparison among the different pumping entities and recharge to GR.The GR pumpage closely followed to the recharge to GR.It is because of the dependence of pumpage to recharge and it this will increase in the future as the demand increases.At the end of the simulation period i.e., at 2150, the demand of PB (combined WP and GR) will be about 10 billion gallons (1% growth).Figure 8 shows the comparison among the different pumping entities and recharge to GR.The GR pumpage closely followed to the recharge to GR.It is because of the dependence of pumpage to recharge and it this will increase in the future as the demand increases.At the end of the simulation period i.e., at 2150, the demand of PB (combined WP and GR) will be about 10 billion gallons (1% growth).

Conclusions
A system dynamics approach facilitated to conduct analysis of the PB watershed and groundwater aquifers with various scenarios.This analysis assisted to estimate WP recharge from water balance and GR recharge via indirect estimation of drawdown.Even though, there were various uncertainties in the parameters (especially in storativity), this study helped to understand some critical future scenarios of PB aquifers.Our study demonstrated that about 90% of water extracted from GR may have been contributed from recharge and GR gets about 0-2 cm/year from recharge (until 2005).The estimated recharge to GR is in the reasonable range while compared to the published values.With this modeling approach, even in an optimistic scenario with larger storativity and larger GR surface area, the historical and current water level depletion may have not been sustained without recharge.Simulation results indicated that the GR will face consistent pressure in the coming years because of the increased growth and probably observe an abrupt decline once recharge threshold is reached.Because of the implied assumptions and parameter values ( with low storativity value), even if the simulation was carried out for longer period i.e., 150 years, the results showed that GR aquifer can face significant drawdown after 2050.As per the previous studies, this study also confirmed that the GR aquifer is at high risk and it is being exploited in an unsustainable way.However, equally, these results should be taken cautiously because various assumptions and constraints had been implied in the model.Our attempt here is rather utilizing the available data than validate the data.The major limitations of this study can be considered as how these estimated recharge values compare to the age dating of groundwater and applicability of the storativity values.

Figure 1 .
Figure 1.Location map of Palouse Basin watershed and surface water and groundwater regions.

Figure 1 .
Figure 1.Location map of Palouse Basin watershed and surface water and groundwater regions.

Figure 2 .
Figure 2. (a) Conceptual schematic of the Palouse Basin Watershed, water-balance based on the variable surface areas and pumping entities (left) (b) Water balance based on the surface area identical to the Palouse Basin (right).

Figure 2 .
Figure 2. (a) Conceptual schematic of the Palouse Basin Watershed, water-balance based on the variable surface areas and pumping entities (left) (b) Water balance based on the surface area identical to the Palouse Basin (right).

Figure 3 .
Figure 3. (a) Palouse Basin groundwater causal loop diagram in a simplest form; (b) Major components of System Dynamics model of the Palouse Basin watershed.

Figure 3 .
Figure 3. (a) Palouse Basin groundwater causal loop diagram in a simplest form; (b) Major components of System Dynamics model of the Palouse Basin watershed.

Figure 4 .
Figure 4. Simulation results of volume of water in (a) Wanapum and Grand Ronde (billion gallons-storativity ´0.0001)(b) curve fitting of drawdown (black for 1% growth, red for 2% growth and purple for some hypothetical midway between 1% and 2%-storativity ´0.001) (c) obseved long-term GR hydrograph at Washington State University (WSU) GR Wells (above mean seal level (a.m.s.l.)).

Figure 5 .
Figure 5. Simulation results of estimated recharge to Grand Ronde and applied drawdown values (cm).

Hydrology 2016, 3 , 13 14 of 19 Figure 6 .
Figure 6.Simulation result of the volume of water in Wanapum and Grande in a billion gallons (1% population growth) assuming no recharge to Grand Ronde.

Figure 7 .
Figure 7. Simulation result of percent contribution of WP to GR as recharge.

Figure 6 . 19 Figure 6 .
Figure 6.Simulation result of the volume of water in Wanapum and Grande in a billion gallons (1% population growth) assuming no recharge to Grand Ronde.

Figure 7 .
Figure 7. Simulation result of percent contribution of WP to GR as recharge.

Figure 7 .
Figure 7. Simulation result of percent contribution of WP to GR as recharge.

Figure 7 .
Figure 7. Simulation result of percent contribution of WP to GR as recharge.

Table 1 .
Period of Availability of Daily Discharge of USGS Gauging Stations of Palouse Basin.

Table 2 .
Mean Areal Recharge of Palouse Basin Sub-Watersheds.

Table 4 .
Initial Volume of Groundwater in Aquifers.