A STOCHASTIC MULTI – SCALE MODEL OF STREAM – GROUNDWATER INTERACTION IN STRONGLY HETEROGENEOUS POROUS MEDIUM AND ITS APPLICATION IN SOUTHERN BRANCH COUNTY, MICHIGAN By Xinyu Ye A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Civil Engineering – Master of Science 2014 ABSTRACT AN INTEGRATED GROUNDWATER RECHARGE AND FLOW MODEL TO PRE A STOCHASTIC MULTI-SCALE MODEL OF STREAM – GROUNDWATER INTERACTION IN STRONGLY HETEROGENEOUS POROUS MEDIUM AND ITS APPLICATION IN SOUTHERN BRANCH COUNTY, MICHIGAN By Xinyu Ye In this paper, stream depletion is assessed by the approach of multi-scale geostatistics in stressed watershed, southern Branch County, Michigan. The watershed is currently under large water demand and representative of the general failure to pass the online Water Withdrawal Assessment Tool. Due to the heterogeneity of porous medium and the high variability of hydrogeological parameters and scale, there is a deviation between field observations and simulated groundwater flow in those areas. The approach of multi-scale geostatistics model based on detailed lithological data and its application in numerical groundwater simulation can be used in stream depletion assessment.     TABLE OF CONTENTS LIST OF TABLES ........................................................................................................................iv LIST OF FIGURES ....................................................................................................................... v 1. INTRODUCTION ..................................................................................................................... 1 Theory of stream depletion ....................................................................................................... 2 Streams and groundwater dependent ecosystem .................................................................... 2 Objectives ................................................................................................................................... 3 2. METHODS ................................................................................................................................. 4 Transition Probability Geostatistics ......................................................................................... 4 Markov Chains Analyses ......................................................................................................... 5 Definition of transition Rates ............................................................................................... 6 Multi-scale Flow Model ............................................................................................................. 7 3. APPLICATION IN SOUTHERN BRANCH COUNTY, MI ............................................... 10 Geology Setting ........................................................................................................................ 13 Conceptual Model .................................................................................................................... 14 Model Development ............................................................................................................... 14 Multi-scale 3D Glacial Aquifer Mapping ........................................................................... 15 Stochastic Multi-Scale quasi-3D Flow Model .................................................................... 18 Data for Multi-scale flow model......................................................................................... 20 Model Calibration .............................................................................................................. 20 4. RESULTS and DICUSSION................................................................................................... 22 Multi-scale 3D Glacial Aquifer Mapping ............................................................................... 22 Regional Aquifer Model -Southern Branch County .............................................................. 22 Local Aquifer Model –Stressed Watershed and Site Scale Model ......................................... 25 Sensitivity Analysis ................................................................................................................ 32 Mean Length Ratio of lateral to vertical ............................................................................ 32 Stochastic Sensitivity Analysis in stability of transition probability geostatistics model ... 34 Stochastic Multi-Scale Quasi -3D Flow Model ...................................................................... 36 Quasi - 3D Regional flow model and local flow model ........................................................ 36 Model Calibration .............................................................................................................. 37 Quasi - 3D site scale flow model ........................................................................................... 37 5.CONCLUSIONS.......................................................................................................................47 REFERENCES.............................................................................................................................48 iii     LIST OF TABLES Table3. 1 Materials setting and conductivity value.....................................................................17 Table3. 2 Details of the multi-scale aquifer model.......................................................................19 Table3. 3 Details of multi-scale flow systems................................................................................19 Table3. 4 Details of well information. .......................................................................................... 20 Table3. 5 Stream leakance and depth in stream order. .................................................................. 20 Table4. 1 Southern Branch County (SBC) characterizations from data........................................23 Table4. 2 Stressed Watershed in SBC, MI characterizations from data. ...................................... 26 Table4. 3. 1 Stressed Watershed in SBC, MI characterizations from data. ................................... 32 Table4. 3. 2 Four different mean length ratios of lateral to vertical. ............................................ 32 Table4.4. 1 Mean value in each Monitoring Wells (MW) in five layers. ....................................... 45 Table4.4. 2 Standard devation in each Monitoring Wells (MW) in each layer ............................ 45 iv   LIST OF FIGURES Figure2. 1 Flow chart of Multi-scale Flow Model. ........................................................................ 8 Figure3. 1 Multi-scale model in regional, local and site scale. .................................................... 11 Figure3. 2 Multi-scale representation of topography for southern Branch County and stressed watershed. ...................................................................................................................................... 12 Figure3. 3 Site scale conceptual model for first - order stream....................................................13 Figure3. 4 Multi-scale lithological data for southern Branch County and Stressed Watershed Area (GWIM 2006) ........................................................................................................................ 16 Figure3. 5 Location of the site scale model and the cross-section of the borehole date. ............. 17 Figure4. 1 Matrix of vertical- direction transition probabilities for Southern Branch County site data. ............................................................................................................................................... 23 Figure4. 2 3D glacial aquifer mapping using the Wellogic data. Resolution in grid NX=285, NY=160, NZ=50; Cell Length X=150m, Y=150m and Z=2.4m. .................................................. 24 Figure4. 3 Cross-sectional view in southern Branch County from West- East and North-South view, unit in meter. ......................................................................................................................... 25 Figure4. 4 Matrix of vertical- direction transition probabilities for stressed watershed site data. ....................................................................................................................................................... 26 Figure4. 5 3D glacial aquifer mapping of stressed watershed using the Wellogic data. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m............................ 27 Figure4. 6 Cross-section in stressed watershed from West-East view, unit in meter. ................... 27 Figure4. 7 Location of the site scale model with first-order stream and the 3D geology mapping using Wellogical data. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m....................................................................................................................... 28 Figure4. 8. 1 3D glacial aquifer mapping of site scale using the Wellogic data and cross sectional view in X direction. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m....................................................................................................................... 29 Figure4. 8. 2 3D glacial aquifer mapping of site scale using the Wellogic data and cross sectional view in X direction. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m....................................................................................................................... 29 v   Figure4. 8. 3 3D glacial aquifer mapping of site scale using the Wellogic data and cross sectional view in X direction. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m....................................................................................................................... 30 Figure4. 9. 1 3D glacial aquifer mapping of site scale using the Wellogic data and cross sectional view in Y direction. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m....................................................................................................................... 30 Figure4. 9. 2 3D glacial aquifer mapping of site scale using the Wellogic data and cross sectional view in Y direction. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m....................................................................................................................... 31 Figure4. 9. 3 3D glacial aquifer mapping of site scale using the Wellogic data and cross sectional view in Y direction. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m....................................................................................................................... 31 Figure4. 10. 1 3D aquifer mapping in case (a). ............................................................................ 33 Figure4. 10. 2 3D aquifer mapping in case (b). ............................................................................ 33 Figure4. 10. 3 3D aquifer mapping in case (c). ............................................................................ 33 Figure4. 10. 4 3D aquifer mapping in case (d). ............................................................................ 34 Figure 4.11 Stochastic analysis in stability of transition probability geostatistics model from ten points. ............................................................................................................................................ 35 Figure4. 12 Head contour and cross-sectional view of southern Branch County and stressed watershed…....................................................................................................................................35 Figure4. 13 Calibration charts for hydraulic heads in southern Branch County.........................36 Figure4. 14 Head contour and cross-sectional view of site scale. The purple color represents the surface drain area. ................................................................................................................... 38 Figure4. 15. 1 Head contour and the conductivity in first computational layer inactive pumping. ....................................................................................................................................................... 38 Figure4. 15. 2 Head contour and the conductivity in second computational layer inactive pumping. ........................................................................................................................................ 39 Figure4. 15. 3 Head contour and the conductivity in third computational layer inactive pumping ....................................................................................................................................................... 39 Figure4. 15. 4 Head contour and the conductivity in fourth computational layer inactive pumping. ........................................................................................................................................ 40 vi   Figure4. 15. 5 Head contour and the conductivity in fifth computational layer inactive pumping. ....................................................................................................................................................... 40 Figure4. 16. 1 Water budget in site scale model inactive pumping. ............................................. 41 Figure4. 16. 2 Water budget in site scale model active pumping. ................................................ 41 Figure4. 17 Location of Monitoring Well (MW) in site scale model.............................................41 Figure4. 18. 1 Condctivity value of Monitoring Wells (MW) in first computional layer..............43 Figure4. 18. 2 Condctivity value of Monitoring Wells (MW) in second computional layer. ........ 43 Figure4. 18. 3 Condctivity value of Monitoring Wells (MW) in third computional layer............. 44 Figure4. 18. 4 Condctivity value of Monitoring Wells (MW) in fourth computional layer. ......... 44 Figure4. 18. 5 Conductivity value of Monitoring Wells (MW) in fifth computational layer. ........ 45 Figure4. 19 Water depletion in first order stream and protected area in each ten realizations....46   vii   1. INTRODUCTION In order to protect the surface-water and groundwater dependent ecosystems, Michigan recently established a new water use law that limits large quantity groundwater withdrawals as part of a concerted effort. The estimated degree of stream depletions by pumping is a significant consideration for ecosystem assessments. Effective implementation of this estimation relies on the ability to quantify stream-aquifer interaction and to estimate pumping induced water depletions in protected area. This, however, proves to be a difficult task because the complexity of surface-water and groundwater systems, glacial sediments are inherently strongly variable. In this paper, a stochastic, multi-scale groundwater model can be used to simulate stream depletion dynamics in strongly heterogeneous media. In particular, we develop a regional and a local and site scale model that can be combined to assess site-specifically the impacts of large quantity withdrawals taking into account strongly aquifer heterogeneity. In addition, transition probability geostatistics is used to create three dimensionally stochastic representation of aquifer geology based on a detailed statewide lithological database. The approach allows simulating geology into the development of geology varibility using iterative Markov chain models (Carle 1999). The output results can be a visual representation on a 3D grid. These material sets can be used for hydrology modeling to simulate site-specifically stream aquifer interaction. The flow model system can be used to simulate quasi – 3D steady flow dynamics in the glacial drift aquifer. 1     Theory of stream depletion Many analytic models have been developed based on the understanding of stream and aquifer systems (Theis, 1941;Glover and Balmer 1954; Hantush 1965; Hunt 1999; Butler et al. 2001). The definition of pumping-induced stream depletion included both induced infiltration of stream water into the aquifer and capture of aquifer discharge to the stream (Theis, 1941; Bredehoeft, 1997; Sophocleous, 1997). Stream depletion rates can be estimated in two-dimensional mathematical models in a fully penetrating stream and pumping well (Theis 1941). In addition, stream aquifer systems in a leaky aquifer system was been discussed (Butler Jr.; Zhan; Zlotnik 2006). The assumptions in most of these theoretical approaches are a uniform aquifer thickness and a continuous streambed of uniform hydraulic conductivity (Hantush 1965). However, the assumptions of horizontal groundwater flow and uniform aquifer limit the applicability in natural systems (Kollet and Zlotnik 2003). Several numerical studies evaluated the parameters’ sensitivity and assessed the impact of various assumptions (Conrad and Beljin, 1996; Butler et al., 2001). Streams and groundwater dependent ecosystem Streams and rivers are open ecosystems which providing nutrients, energy and water (Karr and Dudley 1981). Rivers and streams also play a remarkable function in their physicochemical and biological features (Wang and Seelbach 2008). Many aquatic organisms, especially fishes, depend on the stream flow for their life cycles (Karr and Dudley 1981). Thus it is important to maintain and manage the water depletion rates in order to protect the fish habitat (Murphy and Koski 1989). In addition, stream are the temperature, flow regime, channel morphology and food availability changes can be altered by groundwater pumping (J.F. Wright and A.D. Berrie 1987). 2   Surface-water and groundwater dependent ecosystems are considered as an entire system, one of the important components of the natural environment which support a variety of communities including plants, animals and other organisms due to the significant influence of groundwater on surface water. Objectives The objective of this study is to select a site, which can help to demonstrate the groundwater and surface water interaction in highly heterogeneous area. In particular, a multi-scale geostatistics model based on detailed lithological data is used to simulate three-dimensionally the glacial geology. The robustness of the model will be assessed with sensitivity analysis. In addition, a series of groundwater models, a regional scale model, a local scale model, and one site scale model will be developed to assess site-specifically the impacts of pumping wells in southern Branch County, Michigan. Systematic hydrology models at the regional, local and site scale allows for simulations of stream depletion in targeted streams with the surrounding protected area, using integrated water budget analysis. The specific objectives of this research are: 1. To construct a series of 3D aquifer models for southern Branch County and assess the performance of geostatistics model in sensitivity analysis. 2. To build a hierarchical flow model system to simulate flow dynamics at regional, local and site scale, then quantify the water depletions of target streams and surrounding protected areas using a site scale zone budget. 3   2. METHODS In order to assess the impact of water depletions of targeted streams and surrounding protected areas, the approach of multi-scale geostatistics model based on detailed lithological data are used to simulate 3D glacial geology. The application of multi-scale hydrology modeling integrated with zone-based water budget analysis can be used in stream depletion assessment. In the following part, the transition probability approach and multi-scale flow model approach will be discussed. Transition Probability Geostatistics The transition probability approach is a modified conditioned simulation method from indicator Kriging. This approach simulates the geology using iterative of Markov chain models. The advantages of transition probability geostatistics compared with traditional indicator Kriging methods are that they are conditioned and interpreted geology information and relationships to simulate juxtapositional tendencies (Walker 2002). The T-PROGS package was developed by Carle (1999) to operate using transition probabilities. In geostastic simulation, the time lag is replaced with a distance (Carle, 1999). The transition probability 𝑡!"   ℎ! is expressed in the following, 𝑡!"   ℎ! = 𝑃𝑟  {(𝑐𝑎𝑡𝑒𝑔𝑜𝑔𝑟𝑦  𝑘  𝑜𝑐𝑐𝑢𝑟𝑠  𝑎𝑡  𝑋 + ℎ! )| 𝑐𝑎𝑡𝑒𝑔𝑜𝑔𝑟𝑦  𝑗  𝑜𝑐𝑐𝑢𝑟𝑠  𝑎𝑡  𝑋   } Where ℎ! is a positive lag separation in the direction 𝜑. 4   Markov Chains Analyses The transition probability approach simulates the geology using iterative of Markov chain models (Carle, 1999). The definition of 1-D transition probability matric 𝑇(ℎ! ) is: 𝑇 ℎ! = exp  (𝑅!   ℎ! ) Where 𝑇 ℎ! 𝑡!! (ℎ! ) ⋮ = 𝑡!! (ℎ! ) ⋯ ⋱ ⋯ 𝑡!! (ℎ! ) ⋮ 𝑡!! (ℎ! ) Each matrix 𝑡!"   ℎ! is defined by 𝑡!"   ℎ! = 𝑃𝑟  {(𝑐𝑎𝑡𝑒𝑔𝑜𝑔𝑟𝑦  𝑘  𝑜𝑐𝑐𝑢𝑟𝑠  𝑎𝑡  𝑋 + ℎ! )| 𝑐𝑎𝑡𝑒𝑔𝑜𝑔𝑟𝑦  𝑗  𝑜𝑐𝑐𝑢𝑟𝑠  𝑎𝑡  𝑋   } Where ℎ! is a lag time in the direction 𝜑. 𝑅!   is a transition rate matrix, shown in the following, 𝑅!   𝑟!!,! ⋮ = 𝑟!!,! ⋯ ⋱ ⋯ 𝑟!!,! ⋮ 𝑟!!,! The proportions of each formation are followed, ! 𝑝! = 1 !!! The row sums in 𝑇 ℎ! are, ! 𝑡!" ℎ! = 1    ⍱𝑗 !!! 5   ! 𝑝! 𝑡!"   ℎ! = 𝑝  ⍱  𝑘 !!! Where, 0 ≤   𝑡!"   ℎ! ≤ 1  ⍱𝑗, 𝑘 There are four alternatives methods to generate the trainsition probability approach from Markov chains: definition of transition rates, embedded transition probabilities, embedded transition frequencies, and maximum entropy factors (T-PROGS Version 2.1). In this research, the definition of transition rates is applied here. Thus the following paper will discuss the definition of transition rates. Definition of transition Rates In order to capture the variability of geology unit, Markov chains can be defined in an N times N matrix of transition rates (Carle, 1999). A transition rate is defined as: 𝑟!",! = !"!"   !  ⍱𝑗, 𝑘 !!! Transition rates are calculated in matrix form: 𝑅!   𝑟!!,! ⋮ = 𝑟!!,! ⋯ ⋱ ⋯ 𝑟!!,! ⋮ 𝑟!!,! Thus, the geology varies are described by transition rate in the transition probability in lag time (Carle, 1999). The transition rate matrix are described as: 𝑟!!,! = − 6   1 𝐿!,! Where 𝐿!,! is the average mean length of material k in direction 𝜑. The row sums of 𝑅!   , ! 𝑟!",! = 0    ⍱𝑗 !!! The column sums, ! 𝑝! 𝑟!",! = 0  ⍱𝑘 !!! Multi-scale Flow Model In order to capture the multi-scale flow dynamics and the impact of pumping well, it is necessary to get the multi-scale modeling information in detail. The general approach used in this research is to model the large area to capture the regional scale characteristics, and to refine the model in the local area, modeling more detailed information in site-scale to resolve the flow dynamics. In summary, the all variability, from regional to local to site-scale are solved (Li and his co-workers Liao et al., 2006). The following equations describe the main model 𝐻 at level 𝑙 on Ω! , and its interaction with its “parent” model on Ω!!! at level 𝑙 − 1, and various “daughter” model Ω!!! at level 𝑙 + 1. Parent model provides boundary and initial conditions, reflecting the coupling with the parent at level 𝑙 − 1. 𝑆! 𝜕𝐻! 𝜕 𝜕𝐻!   = 𝐾 𝑥 + 𝑒! 𝜕𝑡 𝜕𝑥 𝜕𝑥 7   BC 𝐻! = 𝐻!!! IC 𝐻! 𝑥, 0 = 𝐻!!! (𝑥, 0) Where, 𝑆 !  is the specific storage, 𝑒 is the source/sink term. BC is the boundary condition. IC is the initial condition. Flow chart of the methods is shown in the Figure 2.1. Figure2. 1 Flow chart of Multi-scale Flow Model. In this research, the stressed watershed which was modeled at the local scale is located in the southwest part in the Branch County, approximately 8000 meters away from the southern Branch County. And for the site scale model, which located in the southern stressed watershed, 9000 meters away from the stressed watershed boundary. The boundary lines for local and site scale 8   model are far away from the site, which will minimize the effect on the results. Thus, the oneway, multi-scale flow model is applied. 9   3. APPLICATION IN SOUTHERN BRANCH COUNTY, MI Branch County is hydrogeologically complex, characterized by strong heterogeneity including structural variation in the vertical direction, highly variable aquifer thickness, very thin drift in localized areas, complex surface topography and presence of ridges and valleys both in the land and bedrock surfaces. Thus Branch County is the general representation of a difficult task because the complexity of surface-water and groundwater systems, glacial sediments are inherently strongly variable. In order to quantify the flow dynamics and the pumping induced water depletion, a series of groundwater models - a regional scale model, a local scale model, and one site scale model to assess site-specifically the impacts of large quantity withdrawals in southern Branch County and stressed watershed, Michigan. Figure 3.1 represents the model set of regional, local and scale model. Complex surface topography of southern Branch County and watershed are shown in Figure 3.2. 10   Figure3. 1 Multi-scale model in regional, local and site scale. 11   Figure3. 2 Multi-scale representation of topography for southern Branch County and stressed watershed. 12   Geology Setting The bedrock geology is dominated by the Michigan Basin, which is an elliptical, intracratonic basin nestled against the southern margin of the Canadian Shield in the Great Lake Region (Robb Gillespie et al.). Branch County is located in the southern most Lower Peninsula of Michigan and is in area of 506 mi2. The St. Joseph watershed drains the county. The topography is variable including flat till plains and rolling hills (Giroux et al. 1966). Based on the February 2005 Wellogic database, about 83 percent of water wells are completed in the glacial deposits, 13 percent are finished in the bedrock and in Branch County. It hard to distinct the last 4 percent wells. Due to their complex depositional history, the glacial deposits in Branch County are lateral and vertical heterogeneous (West john and others 1994). The glacial deposits are up to 500 ft thick in the county. In Michigan, most of the glacial aquifers are sand and gravel (Westjohn and et al., 1994). Wells developed in outwash yield from 10 to 40 gal/min, while, wells finished in till plains produce from 0 to 15 gal/min. Wells developed in thick, coarse outwash deposits can provide up to 5,000 gal/min (Giroux et al. 1966). Based on the Public Water Supply database, the estimated transmissivity ranges from about 590 to 38,700 ft2/day for sand and gravel formations in Branch County (Giroux et al. 1966). Beneath the glacial deposits is Coldwater Shale in Branch County. The thickness of the Coldwater Shale ranges from 440 to 1,000 ft. The formation is mostly blue to gray shale, mixed with thin lenses of limestone, dolomite, and sandstone (Giroux and others, 1966). However, 13   water wells are rarely located in the Coldwater Shale in Branch County. The quantities of water are provided small supply from the wells (Apple and Reeves, 2007). Conceptual Model Since the pumping well is located close to the first order stream in the site scale model, the zone based water budget is only applied in there. At regional scale, southern Branch County provides boundary for stressed watershed. The boundary sets as one way head boundary in flow model, since the boundary of stressed watershed is far away from the southern Branch County. For the site scale model, the stressed watershed provides boundary for it. The boundary also sets one way head boundary in flow model. The site scale conceptual model for this system is shown in Figure 3.3. The multi-scale approach which can integrating the geology and topography features at the regional, local and site scale in order to quantify the flow dynamics. Figure3. 3 Site scale conceptual model for first - order stream. Model Development A multi-scale framework is used to develop a hierarchy of interaction of groundwater and surface water in regional, local and site scale model. The regional model is simulated at a coarse scale 14   enables to resolve the large scale characteristics. At the local model and site scale model, the interested area is selected and modeled in a finer grid. Multi-scale aquifer mapping and flow model are discussed in the following. Multi-scale 3D Glacial Aquifer Mapping In this project, a series of 3D models for the glacial land systems in southern Branch County were constructed using the enhanced Wellogic – a statewide database containing detailed, lithologic information of Michigan’s glacial sediments (Figure 3.4). This approach can characterize the dominant geological structure, both in the horizontal and vertical directions. From the cross-sectional view, it is general representation of a difficult task to simulate because the complexity glacial sediments which are in inherently strongly variable. Figure 3.5 shows the location of the site scale model and the cross-section of the borehole date. 15   Figure3. 4 Multi-scale lithological data for southern Branch County and Stressed Watershed Area (GWIM 2006). 16   . Figure3. 5 Location of the site scale model and the cross-section of the borehole date. Specifically, transition probability geostatistics approach is used to simulate three-dimensionally the glacial geology in southern Branch County. Four kinds of materials are obtained [e.g., AQ (aquifer), MAQ (marginal aquifer), PCM (partially confining material), and CM (confining material)] and presented in a visual representation on a 3D grid. Typical hydraulic conductivity values are assigned for the four kinds of materials. The value of hydraulic conductivity (K) for each material is shown in Table 3.1. These output material can be used for groundwater modeling – to simulate site – specifically stream aquifer interaction. 17   Materials ID 1 2 Formation Kx (m/d) Kx/Ky Kx/Kz Aquifer 30.48 1 10 Marginal 1.524 1 10 Aquifer 3 Partially 0.003048 1 10 Confining Material 4 Confining 0.0003048 1 10 Material Table3. 1 Materials setting and conductivity value. A series of 3D models for the glacial aquifer systems in southern Branch County is using the Wellogic – a statewide database containing detailed, 3D lithologic information of Michigan’s glacial sediments – is applied. Detailed information of regional and local aquifer model is listed in the Table 3.2. Length in X, Y and Z DX DY DZ (m) (m) (m) (m) 42750, 24000 and 122 150 150 2.4 Regional Aquifer Model (southern Branch County) Local Aquifer Model 16382,14815 and 80 80 80 (Stressed Watershed) Table3. 2 Details of the multi-scale aquifer model. 1.6 Stochastic Multi-Scale Quasi-3D Flow Model The regional model will simulate groundwater flow dynamics on a relatively coarse grid in the entire southern Branch County. The local scale model in stressed watershed is able to simulate more detailed flow dynamics. The site-scale models integrated with water budget is able to analysis the pumping induced water depletion in study area. The regional, local, and site scale 18   models will be linked as one way. The regional model will provide boundary conditions to the local model, which the local model will provide the boundary conditions to the site scale. Stream depletion in each vulnerable stream is quantified as a function of time through an integrated, zone-based water budget analysis. Details of the modeling system are provided in Table 3.3. Each realization of transition probability geostatistics model is applied in one flow simulation in stressed watershed. Due to the limited time, total number of ten realizations is used in the stochastic flow analysis. In sum, the stochastic multi-scale 3D numerical flow model between surface water and groundwater is employed to model the pumping reduced water depletion. The model system is run in steady state which represents a long-term average condition in flow system. A steady state model can be considered as sufficient to model the flow system, however for site scale model, it had to be assessed. The pumping well in site scale is 70gpm (381m3/d), which is the minimum value in large quantity withdrawal; detailed information is listed in Table 3.4. Regional flow model Local flow model Site scale flow model Water budget The model domain is in southern Branch County in 38520m ×18818m. The grid resolution is approximately 137 m in one computational layer. The regional flow model is calibrated with the steady state water levels from 732 monitoring wells. The whole flow system is integrated with 10 m DEM. The local flow model is in the stressed watershed scale in 17132m × 15939m. The grid resolution is approximately 100 m in 5 computational layers. The local model will provide the boundary conditions for site scale model. Site scale model in 1650m ×1780m will be used in analyzing the impacts of pumping reduced water depletion. The grid resolution is approximately 10 m in 5 computational layers. This model will be used to investigate pumping well induced impact on the adjacent surface water systems and quantify the total seepage fluxes into/out of surface water bodies and how these fluxes change in response to a stream. Table3. 3 Details of multi-scale flow systems. 19   Well Meter Screen Meter Location Interval x 575357 Top 264 y 140366 Bottom 244 Table3. 4 Details of well information. Data for Multi-scale flow model Since topography is critical to the flow model system. Therefore, 10 m DEM data (USGS, 2006) was used in the flow model system in order to capture the scale variability as well as the regional and local scales. Streams and lakes are modeled as two-way, head dependent model, while land surface elevation is modeled as one-way drains. The leakance of streams in different size is listed in the Table 3.5. Also the Michigan State-wide groundwater database (GMIM 2006) provides the lithology data which can be used in transition probability geostatistics model. For southern Branch County, 12274 points data is used in transition probability geostatistics model. About 8431 points data is applied in watershed transition probability geostatistics model. Recharge for the model systems is assigned as raster from the state-wide database (GWIM 2006). Stream Order Leakance Depth (m/d) (m) 1st 1 0.3 2st 2 0.5 3th 5 1 4th 10 1.3 5th 20 1.6 6th 50 2 Table3. 5 Stream leakance and depth in stream order. Model Calibration The calibration process is necessary to enhance the consistency between the numerical model and the observed model (Anderson and Woessner, 1992). In this study, recharge rate, leakances 20   of rivers and streams, drain leakance and hydraulic conductivity (K) can be modified in the calibration process. A trial - and - error technique was used in the calibration process in order to minimize the calibration statistics, which including the mean absolute residual head, and the root mean squared error (RMSE). The mean absolute residual head describes how well the model is performed overall. The RMSE is calculated by taking a square root of the sum of the square of errors for the observations, defined as: ! Mean absolute residual head = !   ! ! |   ℎ!   − ℎ! ! | RMSE= !   ! ! ! |   21   ℎ!   − ℎ! ! | !.! 4. RESULTS and DICUSSION In this part, the multi-scale three-dimensionally aquifer representations are presented in regional and local scale model. Then sensitivity analyses about the input parameter of mean length ratio of lateral to vertical extent of materials and the stochastic sensitivity analysis of the stability of transition probability geostatistics model are discussed. In addition, the results of stochastic multi-scale three-dimensionally flow model integrating with the water budget of target stream and the limitations are also discussed. Multi-scale 3D Glacial Aquifer Mapping Regional Aquifer Model - Southern Branch County The southern portion of Branch County served as the regional model domain. The transition probabilities for each material can be calculated from the transition rate matrix (Carle 1999). −0.102 0.020 0.263 0.055 𝑅! = 0.100 −0.309 0.071 0.114  𝑚!! 0.104 0.023 −0.193 0.062 0.084 0.020 0.025 −0.131 The regional scale site deposits with a large proportion of aquifer. The formations of four kinds of materials in regional scale model are listed in the Table 4.1. 22   Materials ID Formation Proportion Average vertical lens (m) 1 Aquifer (AQ) 0.437 5.166 2 Marginal Aquifer (MAQ) 0.057 2.728 3 Partially Confining Material (PCM) 0.168 4.804 4 Confining Material (CM) 0.336 6.748 Table4. 1 Southern Branch County (SBC) characterizations from wellogic data.   The average vertical lens length for each material can be calculated. The transition probability curves which generated from the Southern Branch County are represented in the Figure 4.1. Figure4. 1 Matrix of vertical- direction transition probabilities for Southern Branch County site data. This Figure 4.1 interprets the transition probability approach. The red curve indicates the Markov chain which calculated from the transition rates, and the blue line indicates the actual measurements from the borehole data. Due to high proportion of confining material (CM), there is a much greater probability of every material transitioning to it. Since the vertical continuity of 23   geology deposition also reflects the lateral continuity. Thus the mean length ratios of lateral to vertical in AQ, MAQ, PCM and CM are 8.7, 8, 10 and 12, respectively. The 3D aquifer representation model of southern Branch County is shown in Figure 4.2. Intuitive threedimensional mapping has provided very useful insights into the groundwater flow systems and significantly enhanced our understanding of complex aquifer-aquifer connections and groundwater and surface water interactions. Each of the material sets will be conditioned to the wellogical data, which shown in cross-sectional views (Figure 4.3). These material sets can be used for groundwater modeling – to simulate site-specifically stream aquifer interaction. Figure4. 2 3D glacial aquifer mapping using the Wellogic data. Resolution in grid NX=285, NY=160, NZ=50; Cell Length X=150m, Y=150m and Z=2.4m.   24   Figure4. 3 Cross-sectional view in southern Branch County from West- East and North-South view, unit in meter. Local Aquifer Model –Stressed Watershed and Site Scale Model One stressed watershed in Southern Branch County (SBC) is selected as the local model. The transition probabilities for each material can be calculated from the transition rate matrix (Carle 1999). The following Table 4.2 shows the local site characterization. −0.102 0.003 0.032 0.068 𝑅! = 0.052 −0.102 0.011 0.038  𝑚!! 0.124 0.003 −0.196 0.069 0.086 0.003 0.019 −0.109 25   Materials ID Formation Proportion Average vertical lens 1 Aquifer 0.439 5.833 2 Marginal Aquifer 0.021 7.043 3 Partially Confining Material 0.162 5.328 4 Confining Material 0.377 7.783 Table4. 2 Stressed Watershed in SBC, MI characterizations from data. The average vertical lens length for each material can be calculated from transition rate. The transition probability curves which generated from the one stressed watershed are represented in the Figure 4.4. Figure4. 4 Matrix of vertical- direction transition probabilities for stressed watershed site data. Similarity, the red curve indicates the Markov chain which calculated from the transition rates, and the blue line indicates the actual measurements from the borehole data. The mean length ratios of lateral to vertical in AQ, MAQ, PCM and CM are 7.6, 7, 10 and 10. One of the ten 26   realizations is introduced as an example to represent the complex geology. Figure 4.5 is the 3D aquifer plot for the stressed watershed. The cross-sectional view from west to east is shown in Figure 4.6. Figure4. 5 3D glacial aquifer mapping of stressed watershed using the Wellogic data. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m. Figure4. 6 Cross-section in stressed watershed from West-East view, unit in meter. Figure 4.7 shows the location of the site scale model and the 3D geology mapping using Wellogical data. The blue ellipse points out the first-order stream which was modeled at the site scale. As the various cross-sectional slides, complexity glacial sediments are in inherently strongly variable in the site scale. Setting the output as four kinds of material in AQ (aquifer), 27   MAQ (marginal aquifer), PCM (partially confining material), and CM (confining material) is enough to capture the strongly variable world. The confining unit above the pumping well may blocks the groundwater flow. From x and y cross section slides (Figure 4.8.1 to Figure 4.8.3 and Figure 4.9.1 to Figure 4.9.3), the materials are relatively continuity. Figure4. 7 Location of the site scale model with first-order stream and the 3D geology mapping using Wellogical data. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m. 28   Figure4. 8. 1 3D glacial aquifer mapping of site scale using the Wellogic data and cross sectional view in X direction. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m. Figure4. 8. 2 3D glacial aquifer mapping of site scale using the Wellogic data and cross sectional view in X direction. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m. 29   Figure4. 8. 3 3D glacial aquifer mapping of site scale using the Wellogic data and cross sectional view in X direction. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m. Figure4. 9. 1 3D glacial aquifer mapping of site scale using the Wellogic data and cross sectional view in Y direction. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m. 30   Figure4. 9. 2 3D glacial aquifer mapping of site scale using the Wellogic data and cross sectional view in Y direction. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m. Figure4. 9. 3 3D glacial aquifer mapping of site scale using the Wellogic data and cross sectional view in Y direction. Resolution in grid NX=216, NY=190, NZ=50; Cell Length X=80m, Y=80m and Z=1.6m. 31   Sensitivity Analysis Mean Length Ratio of lateral to vertical In transition probability geostatistics model, the vertical transition probability rate can be calculated directly from the formations of four kinds of materials. However, the lateral transition probability is hard to calculate due to insufficient borehole data. It is a huge effect on the aquifer mapping results due to the ratio of lateral to vertical. An example using stressed watershed in Branch County is to do the ratio text. Thus, the following Table 4.3.1 and Table 4.3.2 are listed the proportion and mean length of materials in vertical and four different mean length ratios of lateral to vertical in order to present the effects of the ratio. Materials ID Formation Proportion Average vertical lens 1 Aquifer 0.439 5.833 2 Marginal Aquifer 0.021 7.043 3 Partially Confining Material 0.162 5.328 4 Confining Material 0.377 7.783 Table4. 3. 1 Stressed Watershed in SBC, MI characterizations from data. # Ratio in AQ Ratio in MAQ Ratio in PCM Ratio in CM a 2 3 2 1 b 5 7 5 7 c 7 7 10 10 d 20 20 20 20 Table4. 3. 2 Four different mean length ratios of lateral to vertical. 32   The following four three-dimensional aquifer mappings (Figure 4.10.1 to Figure 4.10.4) represent in case (a), case(b), case(c) and case(d). Figure4. 10. 1 3D aquifer mapping in case (a). Figure4. 10. 2 3D aquifer mapping in case (b). Figure4. 10. 3 3D aquifer mapping in case (c). 33   Figure4. 10. 4 3D aquifer mapping in case (d). Based on the simulation results, the ratio of lateral to vertical is too small in case (a), so that the distribution is so random, not continuity. The geology is more realistic when the ratio is increased in case (b) and (c). However, the geology is too continuity if the ratio is too large in case (d). Stochastic Sensitivity Analysis in stability of transition probability geostatistics model In order to analyze the stochastic characteristics of the realizations, eight points are random picked up from 20, 40, 80, 100, 120, 140, 160, 180, 200, 400, and 600 realizations (Figure 4.11). The tendency of stability of transition probability geostatistics model is increasing with the increasing realizations. The realizations should be enough in order to capture the stochastic characteristics of realizations. For future stochastic analysis of flow dynamics, more than 400 realizations should be needed in stochastic flow simulation. Due to the limited time, 10 each realizations will be applied in the flow model. 34   Figure 4.11 Stochastic analysis in stability of transition probability geostatistics model from eight points. 35   Stochastic Multi-Scale Quasi -3D Flow Model Quasi - 3D Regional flow model and local flow model The regional model domain is in southern Branch County. The grid resolution is approximately 137 m in one computational layer. Figure 4.12 represents the head contour and the crosssectional view in southern Branch County. The head contour and the cross-sectional view illustrate highly variable aquifer thickness, very thin drift in localized areas, and complex surface topography in southern Branch County. The regional model will provide the boundary conditions to the local model. The local model is used to simulate 3D steady groundwater flows at local stressed watershed scale, which designed in 100m grid size in 5 computational layers. Figure4. 12 Head contour and cross-sectional view of southern Branch County and stressed watershed. 36   Model Calibration The regional model is calibrated with the steady state water levels from 732 monitoring wells from the USGS stream-flow gaging stations. Figure 4.13 represents the correspondence between the modeling head and observed head. In the regional model calibration, RMSE is 3.1373, which is acceptable. Figure4. 13 Calibration charts for hydraulic heads in southern Branch County. Quasi - 3D site scale flow model The site scale model is used to simulate quasi 3D steady groundwater flows about the first order stressed stream, which designed in 10m grid size in 5 computational layers. The boundary head is driven from regional model. In the site scale model, 10 different realizations of transition probability approach are used in simulated site scale flow model. As an example, the flow model which used in one transition probability realization is introduced in the following. Figure 4.14 represents the head contour in the first computational layer. The head counter varies due to the geology setting in each computational layer. Figure 4.15.1 to Figure 4. 15.5 show the head contour and the conductivity in each computational layer inactive pumping. From the 37   conductivity map, clay is above the pumping well. Since the interval location of well is 264m to 244m in Z direction, the pumping well pumps water from third layer. Figure4. 14 Head contour and cross-sectional view of site scale. The purple color represents the surface drain area. Figure4. 15. 1 Head contour and the conductivity in first computational layer inactive pumping. 38   Figure4. 15. 2 Head contour and the conductivity in second computational layer inactive pumping. Figure4. 15. 3 Head contour and the conductivity in third computational layer inactive pumping 39   Figure4. 15. 4 Head contour and the conductivity in fourth computational layer inactive pumping. Figure4. 15. 5 Head contour and the conductivity in fifth computational layer inactive pumping. Since the most vulnerable area is the cold water stream and the associated drainage area, the water budget of the site scale model is taking into account of the whole protected area including the drainage area and the first order stream. The following Figure 4.16.1 and Figure 4.16.2 are the site scale zone water budget inactive pumping and active pumping in 381 m3/d (70gpm) 40   using one realization of transition probability model. The results indicate that the water depletion in river and protected area is 100 m3/d, in 26% water which provides to pumping well. In addition, the boundary in the stream provides 213 m3/d water, in 55%. The rest water is coming from rain area. The constant head in and out are assigned as prescribed head from parent model. The error value is acceptable. Transient recharge data can result in high resolution and more accurate input parameter and affect the water budget. Since the 55% water is coming from boundary, the two-way multi- scale model can be consider to assess the effects of pumping well in future. Water Zone Budget inactive pumping Boundary In through river Boundary Out through river River losing water River gaining water Recharge Surface Drain Constant Head In from parent model Constant Head Out from parent model Error Flux (m3/d) 3562 4 874 1914 2813 2042 9473 12763 1 Figure4. 16. 1 Water budget in site scale model inactive pumping. Water Zone Budget active pumping Boundary In through river Boundary Out through river River losing water River gaining water Recharge Surface Drain Constant Head In from parent model Constant Head Out from parent model Pumping well Error Flux (m3/d) 3775 4 874 1819 2814 1928 9473 12764 381 1 Figure4. 16. 2 Water budget in site scale model active pumping. 41   In stochastic analysis, five Monitoring Wells (MW) are selected to capture conductivity value surrounding the pumping well (Figure 4.17 ). The Figure 4.18.1 to Figure 4.18. 5 show the conductivity value in differnet computional layer. It is clearly to see that the variability of condctivity value in five Monitoring Wells (MW) is relatively large in first and last computional layers. In second and thrid computional layers, the condctivity value in different realization varies in smaller range compare with the fisrt or last computional layers. The results of water depletion may relatively close due to the formations of geology setting in third layer in each realization are close. The mean value and the standard devation in each Monitoring Wells (MW) in each layer are listed in the Table 4.5.1 and Table 4.5.2. From the mean value and standard devation of the Monitoring Wells, the fluctation of each ten realization is not too much. The stochastic characteriscts will reflect the water depletion results. Figure4. 17 Location of Monitoring Well (MW) in site scale model. 42   Conduc'vity  in  m/d     35   MW1   30   MW2   25   MW3   MW4   20   MW5   15   10   5   0   1   2   3   4   5   6   7   8   9   10   Each  Realiza'on     Figure4. 18. 1 Conductivity value of Monitoring Wells (MW) in first computional layer. Conduc'vity  in  m/d   30   MW1   25   MW2   20   MW3   MW4   15   MW5   10   5   0   1   2   3   4   5   6   7   8   9   10   Each  Realiza'on   Figure4. 18. 2 Conductivity value of Monitoring Wells (MW) in second computional layer. 43   35   系列1   30   Conduc'vity  in  m/d     系列2   系列3   25   系列4   20   系列5   15   10   5   0   1   2   3   4   5   6   7   8   9   10   Each  Realiza'on     Figure4. 18. 3 Conductivity value of Monitoring Wells (MW) in third computional layer. Conduc'vity  in  m/d     40   系列1   35   系列2   30   系列3   25   系列4   20   系列5   15   10   5   0   1   2   3   4   5   6   7   8   9   10   Figure4. 18. 4 Conductivity value of Monitoring Wells (MW) in fourth computional layer. 44   40   Conduc'vity  in  m/d     35   30   25   20   MW1   15   MW2   10   MW3   5   MW4   0   1   2   3   4   5   6   7   8   9   10   MW5   Each  Realiza'on     Figure4. 18. 5 Conductivity value of Monitoring Wells (MW) in fifth computational layer. Mean MW1 MW2 MW3 MW4 MW5 Layer 1 16.60 11.16 8.85 7.37 13.84 Layer 2 19.09 20.21 21.20 21.29 22.11 Layer 3 19.70 22.16 22.74 24.12 23.9 Layer 4 21.39 23.28 20.02 17.94 22.03 Layer 5 25.15 16.76 17.90 16.80 22.23 Table4.4. 1 Mean conductivity value in each Monitoring Wells (MW) in five layers. SD MW1 MW2 MW3 MW4 MW5 Layer 1 5.1 7.2 17.4 7.7 3.1 Layer 2 6.8 5.9 8.1 8.4 3.0 Layer 3 5.1 7.9 6.9 6.4 3.4 Layer 4 7.3 7.6 7.3 7.8 7.6 Layer 5 6.6 8.7 10.4 9.3 3.6 Table4. 4. 2 Standard deviation of hydraulic conductivity in each Monitoring Wells (MW) in each layer In addition, the results of water depletion in first order stream and protected area of other 9 realizations of transition probability model are summarized in the following Figure 4.19. The mean value of water depletion data is 100.6 m3/d, and standard deviation is 9.0. Due to the regular lithology data distribution, the fluctuation of water depletion in first- order stream and 45   surrounded protected area is relatively small. Although the ten realizations of transition probability model are not enough, the results can still show the fluctuation of the stochastic characteristics. In the next step, more than 400 realizations should be applied in flow model. 140   120   114   Flux  m3/d   100   80   102   102   106   5   6   101   104   102   103   92   80   60   40   20   0   1   2   3   4   7   8   9   10   Each    Realiza'on   Figure4. 19 Water depletion in first order stream and protected area in each ten realization. 46   5. CONCLUSIONS In this research project, the interaction between surface water and groundwater in highly heterogeneous area is investigated by the approach of stochastic multi-scale modeling systems. In particular, the impact of groundwater dependent ecosystems due to pumping well induced water depletion in protected area, is analyzed in zone based water budget. Multiscale transition probability geostatistics approach is applied to simulate three-dimensionally the glacial geology in southern Branch County and stressed watershed. In addition, multi-scale groundwater flow systems is used to model flow dynamics. Based on the results, the pumping induced water depletion in protected area varies in different transition probability realizations. For future work, the two-way multi- scale model will be applied to assess the effects of pumping well. About stochastic analysis of flow dynamics, more than 400 realizations should be needed in stochastic simulation. In transition probability approach, more sufficient borehole data is needed to calculate the lateral probability. The pumping reduced water depletion is the most critical issue which affects the whole groundwater dependent ecosystems, especially for the first order stream ecosystems. Since the first order stream is the most fragile ecosystems related to cold-water habitats, it is critical to protect this vulnerable area. The hierarchical, multi-scale modeling system, which demonstrated in this research, provides a reasonable distribution of materials in aquifer medium, improving numerical groundwater modeling in assessing water depletion in streams and surrounded protected area. 47   REFERENCES 48   REFERENCES Anderson and Woessner, 1992. Quality Assurance and Quality Control in the development and application of ground – water models. Beth A. Apple and Howard W. Reeves, 2007. Summary of Hydrogeologic Condition for the State of Michigan. USGS Report 2007-1236. Butler, J.J. Jr, Tsou, M.S., 1999. The StrpStrm Model for calculation of pumping induced drawdown and stream depletion. Kansas Geol. Survey Computer Ser. Rep. 99-1. Butler, J.J., Zlotnik, V.A., Tsou, M.S., 2001. Drawdown and stream depletion produced by pumping in the vicinity of a finite-width stream of shallow penetration. Ground Water 39(5), 651–659. Cardenas, M.B.R., Zlotnik, V.A., 2003. Three-dimensional model of modern channel bend deposits. Water Resour. Res. 39 (6),1141, doi: 10.1029/2002WR001383, 2003. Conrad, L.P., Beljin, M.S., 1996. Evaluation of an induced infiltration model as applied to glacial aquifer systems. Water Resour. Bull., Am. Water Res. Ass. 32 (6), 1209–1219. Deutsch, C.V., and Journel, A. J., 1992, Geostatistical software library and user’s guide: Oxford University Press, New York, 340 p. Glover, R.E., Balmer, C.G., 1954. River depletion resulting from pumping a well near a river. Eos Trans. Am. Geophys. Union 35, 468–470. Giroux, P. R., Stoimenoff, L. E., Nowlin, J. O., and Skinner, E. L., 1966, Water resources of Branch County Michigan: Michigan Geological Survey Division Report W 6, 158 p. GWIM. Groundwater Inventory and Map Project. State of Michigan Public Act 148. Hantush, M.S., 1965. Wells near streams with semipervious beds. J. Geophys. Res. 70, 2829– 2838. Huasheng Liao, Yang Li, Prasanna V Sampath, Shu-guang Li. Hierarchical Modeling of a Groundwater Remediation Capture System in Michigan. Manuscript 2006. Hill, M.C., 1998. Methods and guidelines for effective model calibration. US Geological Survey, Water-Resources Investigation Report 98-4005, pp. 90. Hunt, B., 1999. Unsteady stream depletion from ground water pumping. Ground Water 37 (1), 98–102. Hunt, B., Weir, J., Clausen, B., 2001. A stream depletion field experiment. Ground Water 39 (2), 283–289. 49   Jenkins, C.T., 1968. Techniques for computing rate and volume of stream depletion by wells. Ground Water 6 (2), 37–46. James J. Butler Jr., Xiaoyong Zhan, and Vitaly A. Zlotnik 2006. Pumping-Induced Drawdown and Stream Depletion in a Leaky Aquifer System. Papers in the Earth and Atmospheric Sciences. Paper 275. James R. Karr and Daniel R. Dudley. 1981. Summary of the black creek project. J.F. Wright and A.D. Berrie 1987.Ecological effects of groundwater pumping and a natural drought on the upper reaches of a chalk stream. Regulated Rivers: Research & Management: Volume 1, Issue 2, pages 145–160, April/June 1987. Justin R. Walker 2002 Application of Transition Probability Geostatistics for Indicator Simulations Involving the MODFLOW Model. K. V. Koski, M. L. Murphy and J. R. Sedell. 1989. Influence of forest practices on aquatic production. Lizhu Wang and Paul W. Seelbach. 2008. A river valley segment classification of Michigan streams based on fish and physical attributes. Transactions of the American Fisheries Society: 137:1621-1636. National Elevation Database. United States Geological Survey (USGS) 2006, http://ned.usgs.gov/. Robb Gillespie, William B. Harrison III, and G. Michael Grammer Michigan Geological Repository for Research and Education Western Michigan University. Sophocleous, M., 1997. Managing water resources systems: why safe yield is not sustainable. Ground Water 35 (4), 561. Steven F. Carle, 1999. T-PROGS: Transition Probability Geostatistical Software Version 2.1. Theis, C.V., 1941. The effect of a well on the flow of a nearby stream. Eos Trans., Am. Geophys. Union 22, 734–738. USGS report understanding and managing the effects of groundwater pumping on streamflow.Jan.2013. Westjohn, D. B., Weaver, T. L., and Zacharias, K. F., 1994, Hydrogeology of Pleistocene Glacial Deposits and Jurassic “Red Beds” in the Central Lower Peninsula of Michigan: U.S. Geological Survey Water-Resources Investigations Report 93-4152, 14 p. Zlotnik, V.A., Huang, H., 1999. Effect of shallow penetration and streambed sediments on aquifer response to stream stage fluctuations (analytical model). Ground Water 37 (4), 599–605. 50   Zlotnik, V.A., Huang, H., Butler Jr., J.J., 1999. Evaluation of stream depletion considering finite stream width, shallow penetration, and properties of streambed sediments. In Proceedings of Water 99, Joint Congress, Brisbane, Australia, pp. 221–226. 51