Water balance between surface water and groundwater in the withdrawal process: a case study of the Osceola watershed

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.

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 km² and of the Chippewa Creek watershed is 7.5 km² (MDEQ 2007). Thus, the drainage area for the stream segment on Twin Creek that will likely be impacted is 43.2 km² while the Chippewa Creek is 7.5 km².

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 m³/d from Chippewa Creek. The combination of the total allowable withdrawal from the two watersheds is 2,616.5 m³/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 m³/h, which is well below the combined allowable withdrawal of 109.0 m³/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 m³/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.

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.

For the test case, the value of T is 92.9 m²/d, d is 152.4 m, S is 0.1 and Q_w is 0.016 m³/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 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.

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.

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

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 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 m³/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.

As indicated in Table 3, the reduced flow quantity in Twin Creek is 408.5 m³/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 m³/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 m³/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 m³/d which is less than 1,417.3 m³/d (the allowable withdrawal from Twin Creek), and the flow intercepted by pumping from Chippewa Creek is 147.6 m³/d which is less than 1,199.2 m³/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.

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.

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.

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.

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).