A generalized hierarchical approach with a water balance function is introduced to simulate stream-flow depletion in a complex groundwater system in Osceola County. The groundwater flow system at the site, because of the complex interaction between ambient streams, exhibits a unique multi-scale pattern that proves to be difficult to simulate using standard modelling tools. The hierarchical modelling system was first calibrated to water level measurements collected from monitoring wells. Afterwards, systematic hierarchical simulations and integrated water budget analyses were performed to evaluate the adverse resource impact of the ongoing water withdrawal. The multi-scale process-based results from this generic hierarchical modelling system provided critical storage and flux information that can be used to comprehensively assess the pros and cons of water resource development and management, such as artificial water withdrawal.

INTRODUCTION

Interaction of surface water and groundwater, as well as its impact on water resources planning and management, is an important and difficult issue to deal with. Groundwater and surface water have often been studied separately, as they have different behaviours both in temporal and spatial scales. However, their flow processes are naturally linked and should be treated as an integrated singular resource when it comes to water planning and management (Bourgault et al. 2014; Boyraz & Kazezyılmaz-Alhan 2014). For example, the base flow in streams and rivers is mostly contributed from groundwater, while pumping water from wells near reservoirs or rivers is most likely recharged from surface water (MacDonald et al. 2014). Land-use conditions play an important role in the interaction of surface water and groundwater (Moiwo & Tao 2014; Bhunia et al. 2014). These physical interaction processes also transport contaminants and thus integrate the biogeochemical interchanges between surface water and groundwater components (Szczucińska 2014). Therefore, water balance modelling needs to include both the contributions from groundwater resources and free surface water resources in a watershed. The basic problem of water balance interactions between groundwater and surface water is stream-flow depletion, i.e., water withdrawals from aquifers can cause adverse impacts on nearby streams, such as stream-flow reduction. Many analytical models have been developed for various aquifer and stream conditions (Wilson 1993; Hunt 1999, 2003; Barbaro & Zarriello 2007; Reeves et al. 2009). Analytical models are easy to use, require minimum data input, offer instantaneous solutions, and are usually assumed to provide conservative estimates. However, for complicated hydro-geologic conditions, numerical models are preferred (Spalding & Khaleel 1991; Sophocleous et al. 1995; Conrad & Beljin 1996; Osman & Bruen 2002; Chen & Shu 2002, 2006; Darama 2004). The complication of a multi-scale flow system for simulating interaction of groundwater and surface water in numerical models is that the groundwater system is mostly in regional scales, while the flow interaction between surface water and groundwater is typically in local scales. Several methods can, in principle, be used to model the multi-scale groundwater and surface flow system. For local grid refinement relatively large size cells in a numerical model grid are subdivided or refined with progressively smaller spatial dimensions in the areas of interest (Gable et al. 1996). An alternative approach to represent detailed well dynamics in a regional model is to use a local analytical correction within the cell containing the well based on the steady-state Thiem equation (Anderson & Woessner 1992). In the approach of local nested numerical correction, grid-dependent information from the regional model is used to construct a separate model with finer grid spacing to obtain more information around the area of interest (Mehl & Hill 2006). The coupled GFLOW-MODFLOW model replaces the upper MODFLOW layer or layers, in which the surface water and groundwater interactions occur, by an analytic element model (GMODFLOW) that does not employ a model grid; instead, it represents wells and surface water directly by the use of point-sinks and line-sinks. Thus it is applicable to relatively large areas, in many cases to the entire model domain, and thus forms an attractive alternative to local grid refinement or inset models (Haitjema et al. 2010). DIVAST, a two-dimensional hydrodynamic and water quality numerical model, was extended to model three-dimensional (3D) groundwater, which integrated the various flow interactions with one model by switching between the shallow water equations and the porous media equations (Sparks et al. 2013). A coupled groundwater and surface flow model was developed for predicting shallow water flows with particular emphasis on its ability to deal with moving boundary problem. A hydrostatic pressure distribution is assumed to enable the water surface gradient to be used as the driving force for the groundwater flow component of the model (Yuan et al. 2008). Huang et al. (2011, 2014) developed a mathematical model for describing 3D groundwater flow induced by a fully penetrating vertical well in aquifers between two parallel streams, in which a general equation is adopted to represent the top boundary condition which is applicable to either a confined, unconfined or leaky aquifer.

The groundwater flow in an area is usually in multi-scale and is characterized by complex interplay of significant variability across disparate length scales. These include variations at ‘well scale’, ‘site scale’ and ‘regional scale’. For large-scale, regional groundwater models with local grid refinement results in a more accurate estimation of hydraulic head or drawdown at the well scale. However, there may be a significant increase in the number of nodes, and the solving process can become very expensive in terms of time and computer resources required. For groundwater models of site scale and well scale, large amounts of detailed field data are needed for determining boundary condition, which are both time-consuming and expensive to obtain. The hierarchical method dispatches the problem of regional scale into many problems of site scale and well scale, and makes it possible to reduce very complex groundwater problems into a number of small, less complex problems that can be solved individually. Li et al. (2006, 2008, 2009) applied the hierarchical approach and water balance function of interactive groundwater (IGW) model to some case studies; however, the focus of the study was not on water budget between surface water and groundwater. Recently, a new module for the interaction between rivers and groundwater was developed and integrated into water balance function of IGW model. However, the effectiveness of this module still needs to be testified. In this study, we employ the hierarchical method, which has been integrated into the 3D IGW model, and is widely used in Michigan State (Liao et al. 2010).

The objective of this study is first to test the 3D IGW model with the new water balance module for estimating the interaction between groundwater and surface water which refers mainly to rivers in large complicated groundwater systems, and then to use this model for an actual aquatic environmental issue in small streams to assess if there is any adverse impact on two trout creeks by a proposed groundwater development project – withdrawal water from pumping wells in the Osceola watershed in Michigan, USA.

MATERIALS AND METHODS

Study area

The study site is located in Osceola Township, which includes Muskegon River, Chippewa Creek and Twin Creek (Figure 1). A company will use groundwater as a source for bottled water, which is a 100% consumptive use with all of the withdrawn water permanently removed from the local hydrologic system. The withdrawal is to be processed via a production well, designated PW-101, and with a maximum proposed pumping rate of 817.7 m3/d. The location of the proposed water withdrawal is the Nestle White-Cedar-Osceola Site (W-C-O-Site) located a quarter mile west of 100th Avenue based on the documents of Michigan State Department of Environmental Quality (MDEQ 2007).

Figure 1

Site map of the study area.

Figure 1

Site map of the study area.

The W-C-O-Site is in the upper reaches of the Chippewa Creek watershed. It is also very close to the Twin Creek watershed. The vicinity of the withdrawal can be characterized as groundwater discharging areas which supply base flow for the two creeks and the Muskegon River. The water well, PW-101, is screened approximately 18.3–42.7 m below the water table (MDEQ 2007). At this depth, the drawdown cone formed by pumping PW-101 will intercept groundwater that naturally discharges further downstream in the two creeks as well as the Muskegon River.

An important component of the adverse resource impact is the drainage area of the affected stream segments. The area of the Twin Creek watershed is 43.2 km2 and of the Chippewa Creek watershed is 7.5 km2 (MDEQ 2007). Thus, the drainage area for the stream segment on Twin Creek that will likely be impacted is 43.2 km2 while the Chippewa Creek is 7.5 km2.

Here an ‘adverse resource impact’ is defined as ‘decreasing the flow of a stream by part of the index flow such that the stream's ability to support characteristic fish populations is functionally impaired’, and the ‘index flow’ is ‘the 50% exceedance flow for the lowest flow month of the flow regime’. The index flow represents a summer low flow condition; it is the median flow for the lowest flow month, and in this case it is in August.

Twin Creek and Chippewa Creek are designated by the Department of Natural Resources (DNR) as trout streams. Historic DNR sampling in these trout streams indicated the presence of numerous brook and brown trout, mottled sculpin, creek chubs and blacknose dace. Few other fish species were found in the samples, which is more evidence of these streams as a trout assemblage. Some supplemental data and analyses also showed these creeks having base flow yields and summer temperatures characteristic of trout streams, and that trout, sculpins and dace were the dominant fish species present.

To determine whether a proposed withdrawal will result in an adverse resource impact, MDEQ must determine what portion of index flow can be withdrawn without functionally impairing the stream's ability to support the characteristic fish populations. According to the DNR's analysis, the maximum portion of the index flow that can be withdrawn at these locations is 13%. More water naturally flows in the stream during other times of the year, except in August when only a small portion, based on summer low flow, can be withdrawn. This protects the stream-flow during the critical low flow period as well as preserves seasonal and annual flow variation, which are important to protect the aquatic habitats for the fish populations. Thus, it is approximately estimated that the allowable withdrawal from Twin Creek is 1,417.3 and 1,199.2 m3/d from Chippewa Creek. The combination of the total allowable withdrawal from the two watersheds is 2,616.5 m3/d.

The proposed withdrawal would intercept some of the groundwater discharging into Twin Creek and Chippewa Creek. Therefore, the influences of the proposed withdrawal are evaluated against the combined allowable withdrawal from both the creeks. The proposed maximum withdrawal is 34.1 m3/h, which is well below the combined allowable withdrawal of 109.0 m3/h. Further, the proposed withdrawal is less than the allowable withdrawal from the watershed of either creek alone. Based on the above information, the MDEQ finds that the proposed withdrawal of 34.1 m3/h from the described W-C-O-Site is not likely to cause an adverse impact on either Twin Creek or Chippewa Creek (MDEQ 2007).

It is very interesting to know how much water will really be withdrawn from each creek. Thus, a 3D groundwater model with water balance functionality is needed to gain the correct water distribution. In this study, the 3D groundwater model (IGW) was introduced to simulate the process of groundwater flow in the watersheds, and then the effect of the withdrawal by pumping well was evaluated based on the predictions of the withdrawal from the two creeks and the Muskegon River. The basic data, including hydrogeology condition of the aquifer and hydrologic condition of creeks and rivers for the simulation, are from the geographic information system (GIS) database of MDEQ in Michigan State. Other parameters, such as hydraulic conductivities, recharges and river leakages, which are used in the simulation, are calibrated from the site pumping tests provided by MDEQ.

METHODS

The 3D model for groundwater simulation (IGW) (Li et al. 2008) is selected to study water balance in a complicated groundwater system. Li et al. (2009) developed a new method of implementing nested grid modelling in the groundwater model. In particular, this method took advantage of two emerging computing strategies – ‘dynamic fusion’ and ‘interactive steering’ – to develop a dynamically integrated, human-centred ‘hierarchical patch dynamics paradigm (HPDP)’ for generalized hierarchical groundwater modelling. The two modern computing strategies make it possible to employ arbitrary numbers of nested sub-models (patch models) for solving much more complex groundwater problems. Instead of treating the simulation of groundwater flow and solute transport separately, it models both processes concurrently (sequential within a time step). Instead of treating the various scales of patch dynamic modelling (regional, sub-regional, local, site, hot spots) as different phases in a long sequential batch-process, the multi-scaled processes are coupled and modelled simultaneously.

Mathematical model

In hierarchical modelling, the governing equation in each model is the same as the general groundwater governing equation. More specifically, assuming that there are L nested model levels and P(l) patches in the lth model level in a multi-scale modelling system denoted as Mpl [P = 1, P(l), l = 1, L], the governing equations in this system, then, are as follows: 
formula
1
where is the aquifer specific storage coefficient ; H is piezometric head ; t is time [T]; ∇ is the gradient operator ; K is the saturated hydraulic conductivity tensor [L/T]; q represents source (positive) or sink (negative) terms including pumping wells, streams, lakes, drains, etc.; Γ1 is the computational domain boundary on which head distribution, f, is known; Γ2 is the boundary where a known flux, g, is specified; is the spatial vector and is head distribution over the whole computational domain at initial time level t = 0.

Hierarchical patch

The HPDP was originally developed by Wu and his colleague to address complexities and scale interactions in the context of landscape ecology and ecological modelling (Wu & Loucks 1995, 1998). This approach makes a complex system represented hierarchically by breaking it ‘vertically’ into many levels and ‘horizontally’ into many patches (Figure 2, left). The dynamics in different patches are related efficiently to each other through a hierarchical patch network (Figure 2, right). The patches at higher levels are used to simulate larger-scale dynamics whereas nested patches at lower levels are used to simulate smaller-scale processes. The upper level (level l-1) exerts constraints (e.g., as boundary conditions) on the lower level (level l + 1). The conceptual basis for the hierarchical approach emerges from hierarchy theory and a diversity of studies in various disciplines, including management science, biology, ecology and system science. The hierarchy theory has been significantly expanded in the context of evolutionary biology and ecology.

Figure 2

Hierarchical decomposition of a groundwater system (left) recursive subdivision and hierarchical patch network.

Figure 2

Hierarchical decomposition of a groundwater system (left) recursive subdivision and hierarchical patch network.

Model test

Hunt (1999, 2003) presented an analytical solution for a well near a river. The conceptual model is shown in Figure 3. Depicted on the left is a river and on the right is an aquifer. The stream edge was modelled as an infinitely long straight line with zero drawdown. The stream was assumed to completely penetrate the homogeneous aquifer. The changing in free surface elevations was assumed to be small enough to allow use of the linear form of the equations. There is a well in the aquifer. Its pumping rate is Qw and its distance from the stream is d. The solution deducted by Hunt is as follows: 
formula
2
where Qs is the rate of stream-flow depletion (length cubed per unit time); Qw is the pumping rate (length cubed per unit time); erfc() is the complementary error function (dimensionless); d is the distance from the well to the stream (length); S is the storage coefficient, or storability, of the aquifer (dimensionless); T is the transmissibility of the aquifer (length squared per unit time); t is the time from the start of pumping.
Figure 3

Concept model of fully penetrating stream with no stream bed resistance.

Figure 3

Concept model of fully penetrating stream with no stream bed resistance.

For the test case, the value of T is 92.9 m2/d, d is 152.4 m, S is 0.1 and Qw is 0.016 m3/s, respectively. In the 3D model the whole computing area is 2.0 km × 2.0 km, the thickness of the aquifer is 92.9 m and k is 1.0 m/day. This renders T equivalent to the value given in the analytical model, that is to say the conditions of the Hunt's model are the same as those of the 3D IGW model.

The results from the 3D model for 100 days of pumping are shown in Figure 4. To test the accuracy of the 3D model, the analytical solution of Hunt's model was also given in the figure. It is noticeable from Figure 4 that the results from the 3D model agree well with those of the analytical solution.

Figure 4

Stream-flow depletions by the 3D IGW model to be tested with the analytical solution under an ideal case in which two solutions should be the same.

Figure 4

Stream-flow depletions by the 3D IGW model to be tested with the analytical solution under an ideal case in which two solutions should be the same.

Modelling calibration

The whole computing area is about 24.1 km × 23.3 km. The recharge in the area is about 20.3–60.8 cm/year. According to the GIS database of MDEQ, the permeability in the horizontal direction around this region is about 1.42–4.29 m/day. The outer boundaries are represented by no-flow boundaries in the regional model, as shown in Figure 5.

Figure 5

Spatial distribution of groundwater level and flow field simulated in the Osceola watershed.

Figure 5

Spatial distribution of groundwater level and flow field simulated in the Osceola watershed.

The groundwater flow model was calibrated under a steady-state condition. The calibration is conducted based on site scale model. The calibration parameters are conductivities, recharges and leakage coefficient of the streams. The calibration targets are static water levels collected at monitoring wells throughout the site and seepage flux from USGS gauge and field tests. Calibration was made by minimizing the sum of squared differences between the simulated and observed heads, and simulated and observed flux. The final values of parameters obtained, such as conductivities, recharge, etc., are presented in Table 1. The simulated results of groundwater tables and flux in the region are shown in Figure 5.

Table 1

Values of parameters selected as model coefficient

 Conductivity (m/day) Recharge (mm/year) Twin Creek leakage coef. Chippwa Creek leakage coef. Muskegon River leakage coef. 
Site 1.42–4.29 47.5–62.2 0.46 0.46 4.65 
 Conductivity (m/day) Recharge (mm/year) Twin Creek leakage coef. Chippwa Creek leakage coef. Muskegon River leakage coef. 
Site 1.42–4.29 47.5–62.2 0.46 0.46 4.65 

As indicated in Figure 5, there are nine monitoring wells where we have collected observed data on the groundwater levels. Therefore, these data were compared with the simulated results, as shown in Table 2, to justify the applicability and accuracy of the model. From Table 2, we can get the average relative error of 0.04%. It is obvious that the simulated results fit with the observed results quite well.

Table 2

Water head comparison between the observed and the simulated at monitoring wells

Well no. X (m) Y (m) Observed head (m) Simulated head (m) Error (m) Relative error (%) 
548,618.9 370,114.0 315.42 315.55 0.13 0.04 
549,429.8 373,436.4 343.57 343.55 −0.02 0.01 
550,405.9 368,459.4 299.21 299.32 0.11 0.04 
558,638.5 376,445.0 332.97 332.98 0.01 0.01 
560,840.5 380,483.5 348.83 348.94 0.11 0.03 
560,963.5 376,066.3 317.62 317.40 −0.22 0.07 
562,958.9 365,924.7 343.41 343.69 0.28 0.08 
563,563.4 373,327.5 312.00 312.01 0.01 0.01 
565,538.2 373,978.5 313.54 313.43 −0.11 0.04 
Well no. X (m) Y (m) Observed head (m) Simulated head (m) Error (m) Relative error (%) 
548,618.9 370,114.0 315.42 315.55 0.13 0.04 
549,429.8 373,436.4 343.57 343.55 −0.02 0.01 
550,405.9 368,459.4 299.21 299.32 0.11 0.04 
558,638.5 376,445.0 332.97 332.98 0.01 0.01 
560,840.5 380,483.5 348.83 348.94 0.11 0.03 
560,963.5 376,066.3 317.62 317.40 −0.22 0.07 
562,958.9 365,924.7 343.41 343.69 0.28 0.08 
563,563.4 373,327.5 312.00 312.01 0.01 0.01 
565,538.2 373,978.5 313.54 313.43 −0.11 0.04 

To further acquire accuracy of the model in regional scale, a relationship, as shown in Figure 6, of groundwater levels obtained by the model and by observation from the monitoring wells of the MDEQ's database was established. As indicated in Figure 6, the root mean square error of the groundwater levels is approximately 7.18 m or 8.7% of the maximum observed head difference across the site. The arithmetic mean error is 1.93 m or 2.3% of the maximum observed head difference.

Figure 6

Comparison of observed and simulated water heads at monitoring wells in the site model.

Figure 6

Comparison of observed and simulated water heads at monitoring wells in the site model.

On the basis of these comparisons, we can find that the simulation model has satisfactory accuracy for the regional groundwater estimation.

RESULTS

To evaluate the effect on the two trout creeks from the proposed withdrawal, we set two scenarios in which one includes the pumping from the well and the other does not. For a closed watershed with no flow boundary, groundwater will not flow into or out from its boundary. Thus in the simulation, it is assumed that no groundwater flows into or out by the boundary in the study region and only infiltration recharges flow into the aquifer. If groundwater that flows out of the watershed only discharges into rivers, the water quantity from infiltration recharges ought to be equal to the total flux which discharges into all creeks and rivers. Therefore, the groundwater model with the hierarchical approach and the water balance function were applied to simulate the groundwater flow field and water budget.

The velocity field and the groundwater level in the second scenario with the pumping well are shown in Figure 7. To gain water budget between the two creeks and the pumping well, a sub-model was introduced to reveal the groundwater flow in more details. Moreover, the sub-model function can display the process of stream-flow depletion for both large areas and small-scale local regions.

Figure 7

Groundwater flows simulated by the main model (left) and the sub-model (right).

Figure 7

Groundwater flows simulated by the main model (left) and the sub-model (right).

Velocity field and groundwater-level contour

The flow field and groundwater-level contour are presented in Figure 7. The simulated results of the main model show complex groundwater flow in the whole study area, while the sub-models provide a relatively steep hydraulic gradient near the pumping well and the two creeks. With regional recharge, the entire flow of groundwater in the aquifer runs from both the north and south to the middle region, and finally into the Muskegon River based on the results of the main model. Relative to the recharge flow, the pumping well intercepts a small quantity of water which is approximately 34.1 m3/h from the ambient area. Thus, the effect of the pumping is concentrated only in the local area and will not change obviously the spatial distributions of the whole groundwater levels in the whole study region.

Through this method, the main model at higher level is used to simulate larger-scale groundwater dynamics and the sub-models at lower levels are used to simulate smaller-scale processes. Therefore, in this way, all of the groundwater flow both in the regional scale and in the local scale can be displayed in the necessary detail.

Water balance in the watershed

Figure 8 shows the process of water balance into and out from the aquifer in the natural conditions, namely without pumping. Positive flow above the lateral axis in Figure 8 indicates that creeks and rivers get a water supply from groundwater, while the negative flow below the lateral axis implies creeks and rivers lost water into the aquifers. Infiltration recharge is distributed by the lost water from Twin Creek, Chippewa Creek, Muskegon River and other rivers. The recharge into the aquifer is in total 3.99 × 105 m3/d. For Twin Creek and Muskegon River, both of them only gain water from the aquifer. The supplying rate from groundwater is 2.62 × 104 and 2.51 × 105 m3/d individually. For Chippewa Creek, parts of its reaches gain water from the aquifer and the gaining flux rate is 9.05 × 103 m3/d, while other reaches supply water into the aquifer and the supplying flux is 2.81 × 103 m3/d. By using the water balance module, some of the other rivers will gain groundwater of 1.25 × 105 m3/d, while the remainder of the other rivers will have lost surface water of approximately 9.74 × 103 m3/d. In summary, the flux rate of their gaining groundwater is 1.15 × 105 m3/d. The error of water budget in the aquifer is 0.7 m3/d and its relative error of recharge is 0.0001%. It is obvious that the error is quite small by applying the water balance function in this simulation. Therefore, the water balance module in the 3D IGW model is capable of dealing with the interaction between groundwater and river systems.

Figure 8

Water balance in the Osceola watershed without pumping condition.

Figure 8

Water balance in the Osceola watershed without pumping condition.

 Figure 9 illustrates the distributions of precipitation recharging through infiltration in the study watershed. Owing to the pumping effect of a proposed pumping rate of 817.7 m3/d, some of the groundwater which originally flows into the two creeks, the Muskegon River and other rivers will be intercepted by the pumping well. Therefore, stream-flow will be affected by the withdrawal of the pumping well. The comparison of water fluxes with and without pumping well is given in Table 3.

Table 3

Flux comparison under the conditions with and without the pumping well

 Twin Creek Chippewa Creek Muskegon River Other Recharge Pumping rate 
Flux (no well, m3/d) 26,578.9 6,239.6 251,306 115,228 399,352. 4 
Flux (well, m3/d) 26,170.4 6,091.9 251,186 115,088 399,352.4 817.7 
Effect (m3/d) 408.5 147.6 120.5 139.9 816.5 (flow into the pumping well) 
Percentage (%) 50.0 18.1 14.7 17.1   
 Twin Creek Chippewa Creek Muskegon River Other Recharge Pumping rate 
Flux (no well, m3/d) 26,578.9 6,239.6 251,306 115,228 399,352. 4 
Flux (well, m3/d) 26,170.4 6,091.9 251,186 115,088 399,352.4 817.7 
Effect (m3/d) 408.5 147.6 120.5 139.9 816.5 (flow into the pumping well) 
Percentage (%) 50.0 18.1 14.7 17.1   
Figure 9

Distribution of infiltration recharge in the Osceola watershed with pumping condition.

Figure 9

Distribution of infiltration recharge in the Osceola watershed with pumping condition.

As indicated in Table 3, the reduced flow quantity in Twin Creek is 408.5 m3/d. This reduction accounts for only 1.5% of the flow in Twin Creek. However, it makes a great contribution to the proposed withdrawal and reaches to approximately 50% of the proposed pumping rate; whereas, the contributions of Chippewa Creek, Muskegon River and other small rivers to the proposed withdrawal are small. Their contributions are 18.1% in Chippewa Creek, 14.7% in Muskegon River and 17.2% in other creeks, respectively.

With a proposed pumping rate of 817.7 m3/d, the simulation results of the water balance module show that the flux flowing into the pumping well from the region is equal to 816.5 m3/s with a relative error of 0.15%, which demonstrates the accuracy of water balance function in this model. It is clear that the flow intercepted by pumping from Twin Creek is 408.5 m3/d which is less than 1,417.3 m3/d (the allowable withdrawal from Twin Creek), and the flow intercepted by pumping from Chippewa Creek is 147.6 m3/d which is less than 1,199.2 m3/d (the allowable withdrawal from Chippewa Creek). Thus, the proposed withdrawal from the described W-C-O-Site is not likely to cause an adverse impact on either Twin Creek or Chippewa Creek.

DISCUSSION

The permeability of the whole study region is different and hydraulic conductivities in the horizontal direction vary from 1.42 to 4.29 m/day. From the hydrogeology GIS database of MDEQ in the Michigan State, hydraulic conductivities were directly imported into the model as initial conditions. After parameters were calibrated in numerical simulation, different zones had different hydraulic conductivities. However, for the analytical model, only one average value of hydraulic conductivity can be used for the prediction of flow depletion in rivers. Moreover, the location of the pumping well is closer to Chippewa Creek than to Twin Creek. If an analytic model is applied, the pumping well will gain more water from Chippewa Creek than Twin Creek. However, from the results of the water budget based on the simulation model, water from Chippewa Creek is less than that from Twin Creek, and only a quarter of water quantities from Twin Creek. Table 4 presents the resulting contrast between the numerical model and the analytical solution. On the basis of the above analysis and due to the complicated conditions in this study, the numerical model will be more preferable and will give a more reasonable result.

Table 4

Flow depletion contrast between the numerical model and the analysis solution

 Twin Creek Chippewa Creek Muskegon River Other rivers 
Numerical model 408.5 m3/d 147.6 m3/d 120.5 m3/d 139.9 m3/d 
Analysis solution 114.1 m3/d 660.9 m3/d 17.0 m3/d 25.6 m3/d 
 Twin Creek Chippewa Creek Muskegon River Other rivers 
Numerical model 408.5 m3/d 147.6 m3/d 120.5 m3/d 139.9 m3/d 
Analysis solution 114.1 m3/d 660.9 m3/d 17.0 m3/d 25.6 m3/d 

The hierarchical method in this model satisfactorily solves the multi-scale groundwater system and the interactional flow system of groundwater and surface water in this study. It dispatches the groundwater system as a parent model and the interactional part as a sub-model. The parent model presents groundwater flow process in regional scale and the sub-model deals with the flow interactions among creeks, rivers and groundwater, including the pumping well. In the process of the simulation, groundwater, exchanging of surface water and groundwater, and impact on groundwater from the pumping are integrated, and detail flows both in large scale and small scale are presented in the parent model and sub-model individually. Their interactional relations with each other are held the same on their overlapping boundary through iterative computation between the parent model and the sub-model. Generally, the multi-scale problem regarding the interaction of subsurface water and surface water is important and difficult to deal with during the simulation, but the proposed hierarchic method in this study can overcome this problem. It deals with the whole computing area of 24.1 km × 23.3 km as a parent model, and the sub-model region of 0.42 km × 0.39 km which represents a region with a complicated groundwater system characterized by the pumping well, two creeks and rivers.

The water balance module in this integrated model system can display the water budget in the whole watershed clearly. Moreover, the water budget between the pumping well and the aquifer, water balance among the aquifer, creeks and rivers are also simulated satisfactorily. Furthermore, on the basis of the simulated results, especially the results of impact evaluation of the proposed withdrawal at the W-C-O-Site pumping well on the trout creeks, the water balance module demonstrates its applicability in complicated river basin conditions.

There are some limitations to this study. First, only the annual water balance was calculated. Generally, the water level in rivers varies within a year because of the uneven distribution of annual precipitation. This is partly because the main purpose of this study is to evaluate the annual effect of proposed groundwater withdrawal on the trout creeks in the Osceola watershed and to verify the applicability of the water balance module in the 3D IGW model. Second, due to lack of sufficient field data, hydraulic conductivities are calibrated only in the large region of the whole catchment, and no detailed values could be determined at the small-scale zone around the pumping well, creeks and rivers. These will affect the result in the numerical model. As well, the leakage coefficients are calibrated only on three points and will have also some offset effect on the results. All of these limitations will cause some errors in the simulation. However, from the results of water balance in the watershed and the pumping well effect, we can find that this simulation model has a satisfactory accuracy.

CONCLUSIONS

In this paper, the generic hierarchical approach and the water balance interaction function between surface water and groundwater are introduced to simulate stream-flow depletion in a complicated groundwater system of the Osceola watershed. The hierarchical modelling system was first calibrated based on the observed water levels collected from the monitoring wells. Then, systematic hierarchical simulations and integrated water budget analyses were performed to evaluate the impact of the ongoing water withdrawal by the pumping well at the W-C-O-Site on the flows in the trout creeks and rivers at annual scale. The specific conclusions are as follows:

  • The hierarchical model can simulate reasonably well the groundwater flow in the complicated groundwater system.

  • The water balance module in this model can display detailed water withdrawal information by the pumping well.

  • The results of water balance modelling demonstrate that the proposed withdrawal of groundwater is not likely to have an adverse resource impact on the trout creeks to maintain their normal ecological functions.

ACKNOWLEDGEMENTS

First, we would like to express our great gratitude to China Scholarship Council for providing financial support for Dr Bing-dong Li's study in the USA. Second, we would like to thank Prof. Shu-guang Li, the director of the laboratory of Excellence for Real time Computing and Multi-scale Modelling, and Prof. Hua-sheng Liao, visiting Professor in the Department of Civil & Environmental Engineering in Michigan State University, for their supervision and valuable suggestions, and for providing us with valuable data. Special thanks are also extended to the two anonymous reviewers. With their valuable and deep insightful comments and suggestions, the qualities of this paper are improved substantially. This study has also been partially supported by the National Natural Science Foundation of China (Grant Nos 51379137, 81200781) and the open funds from State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University (SKLH-OF-1101and SKLH-OF-1212).

REFERENCES

REFERENCES
Anderson
M. P.
Woessner
W. W.
1992
Applied Groundwater Modeling: Simulation of Flow and Advection Transport
.
Academic Press
,
San Diego, CA
.
Barbaro
J. R.
Zarriello
P. J.
2007
A precipitation-runoff model for the Blackstone River Basin
,
Massachusetts and Rhode Island: US Geological Survey Scientific Investigations Report, Reston, VA, 2006–5213
.
Bourgault
M. A.
Larocque
M.
Roy
M.
2014
Simulation of aquifer-peatland-river interactions under climate change
.
Hydrol. Res.
45
(
3
),
425
440
.
Gable
C. W.
Trease
H. E.
Cherry
T. A.
1996
Geological applications of automatic grid generation tools for finite elements applied to porous flow modeling
. In:
Proceedings of the 5th International Conference on Numerical Grid Generation in Computational Fluid Dynamics and Related Fields
,
Engineering Research Center, Mississippi State University Press, MS, USA
.
Haitjema
H. M.
Feinstein
D. T.
Hunt
R. J.
Gusyev
M. A.
2010
A hybrid finite-difference and analytic element groundwater model
.
Ground Water
48
,
538
548
.
Li
S. G.
Liu
Q.
2006
A real-time, computational steering environment for integrated groundwater modeling
.
Ground Water
44
,
758
763
.
Li
S. G.
Liao
H. S.
Abbas
H.
2008
Computational discovery and innovation in groundwater science and engineering, MODFLOW and more 2008
.
International Groundwater Modeling Center, Colorado School of Mines
,
Golden, Colorado, USA
.
Li
S. G.
Liao
H. S.
Afshari
S.
Oztan
M.
Abbas
H.
Mandle
R.
2009
A GIS-enabled hierarchical modeling patch dynamics paradigm for modeling, complex groundwater systems across multiple scales. Invited book chapter
. In:
Modeling of Pollutants in Complex Environmental Systems
,
Vol. 1
(
Hanrahan
G.
, ed.).
ILM Publications
,
St, Albans
,
UK
, pp.
177
216
.
Liao
H. S.
Li
S. G.
Li
Y.
2010
GIS-enabled multiscale modeling of a complex groundwater remediation site in Michigan
. In:
AWRA Spring Specialty Conference. Geographic Information Systems and Water Resources VI
.
Orlando, FL
.
MacDonald
A. M.
Lapworth
D. J.
Hughes
A. G.
Auton
C. A.
Maurice
L.
Finlayson
A.
Gooddy
D. C.
2014
Groundwater, flooding and hydrological functioning in the Findhorn floodplain, Scotland
.
Hydrol. Res.
45
,
755
773
.
Michigan Department of Environmental Quality
2007
ICE Mountain Spring Water Petition for a No Adverse Resource Impact Determination
. .
Reeves
H. W.
Hamilton
D. A.
Seelbach
P. W.
Asher
A. J.
2009
Ground-water-withdrawal component of the Michigan water-withdrawal screening tool: US Geological Survey Scientific Investigations Report 2009–5003
.
Sophocleous
M.
Koussis
A.
Martin
J. L.
Perkins
S. P.
1995
Evaluation of simplified stream–aquifer depletion models for water rights administration
.
Ground Water
33
,
579
588
.
Sparks
T. D.
Bockelmann-Evans
B. N.
Falconer
R. A.
2013
Laboratory validation of an integrated surface water-groundwater model
.
J. Water Resour. Protect.
5
,
377
394
.
Wilson
J. L.
1993
Induced infiltration in aquifers with ambi­ent flow
.
Water Resour. Res.
29
,
3503
3512
.
Wu
J.
Loucks
O. L.
1998
Hierarchical patch dynamics as a framework for scaling
. In:
Proceedings of the Workshop on Scaling and Modeling in Forestry, Applications in Remote Sensing and GIS
,
University of Montreal
,
Montreal
, pp.
64
71
.