Simulating streamflow and water table depth with a coupled hydrological model

Alphonce Chenjerayi GUZHA*, Thomas Byron HARDY

1. Department of Agricultural and Biological Engineering, Southwest Florida Research and Education Center, University of Florida, Immokalee, FL 34142, USA

2. River Systems Institute, Texas State University, San Marcos, Texas 78666, USA

Simulating streamflow and water table depth with a coupled hydrological model

Alphonce Chenjerayi GUZHA*1, Thomas Byron HARDY2

1. Department of Agricultural and Biological Engineering, Southwest Florida Research and Education Center, University of Florida, Immokalee, FL 34142, USA

2. River Systems Institute, Texas State University, San Marcos, Texas 78666, USA

A coupled model integrating MODFLOW and TOPNET with the models interacting through the exchange of recharge and baseflow and river-aquifer interactions was developed and applied to the Big Darby Watershed in Ohio, USA. Calibration and validation results show that there is generally good agreement between measured streamflow and simulated results from the coupled model. At two gauging stations, average goodness of fit (R2), percent bias (PB), and Nash Sutcliffe efficiency (ENS) values of 0.83, 11.15%, and 0.83, respectively, were obtained for simulation of streamflow during calibration, and values of 0.84, 8.75%, and 0.85, respectively, were obtained for validation. The simulated water table depths yielded averageR2values of 0.77 and 0.76 for calibration and validation, respectively. The good match between measured and simulated streamflows and water table depths demonstrates that the model is capable of adequately simulating streamflows and water table depths in the watershed and also capturing the influence of spatial and temporal variation in recharge.

hydrological modeling; model coupling; streamflow; groundwater; TOPNET model; MODFLOW model; Big Darby Watershed

1 Introduction

In recent years, integrated surface-subsurface modeling tools have evolved rapidly and are now being applied in various watershed hydrology studies in different parts of the world. Studies include those by Markstrom et al. (2008), Panday and Huyakorn (2004), Jones et al. (2006), and Werner et al. (2006). The application of coupled models provides evidence of the capacity of these models to produce realistic catchment behavior. However, as stated by Nemeth and Solo-Gabriele (2003), linking groundwater and surface water models is frequently problematic because the models use different sets of governing equations. Additionally, the time scale is usually longer for groundwater modeling than for surface water modeling.

Considerable effort has been expended to characterize the physical, chemical, and biological processes affecting groundwater and surface water resources in river basins. This is because it hasbecome apparent in hydrological studies that processes must be perceived in an integrated way. Many of the impacts of land use changes on surface water systems cannot be evaluated meaningfully without considering the dynamics in subsurface flow systems. As the development of fully integrated model concepts for this purpose is still in its early stages, one means of integration is the coupling of existing models. This leads to other problems, however, because most models have been designed to simulate specific aspects of the water cycle. Coupling of two or more models can also result in inconsistencies because the individual models may describe the same processes in different ways. In various studies, the coupling of surface and subsurface flow models has begun with the setup of relationships between river stages and groundwater storage (Pinder and Sauer 1971). Early attempts at coupling hydrological models included MODBRANCH (Swain and Wexler 1993), which couples the groundwater flow model MODFLOW and the river network program BRANCH. More recently, Smits and Hemker (2004) modeled the interaction of surface water and groundwater flow by linking Duflow to Microflow. Ellingson and Schwartzman (2004) integrated an unsaturated zone flow model and a groundwater model in the regional HSPF Model.

In 2008, the U. S. Geological Survey released another coupled model, GSFLOW (Markstrom et al. 2008), which integrates MODFLOW and the precipitation-runoff modeling system (PRMS). This approach of coupling already existing models will likely reduce costs in development of new integrated model codes. In most cases where model integration has been attempted, it is very costly to build a single predictive model that adequately represents all hydrological processes, and it is therefore important to link models of individual processes. This research was inspired by the need to improve tools for simulating interactions between groundwater and surface water to quantify the effects of human activity and natural phenomena on watershed hydrological responses. The research was carried out using the TOPNET and MODFLOW models with application to the Big Darby Watershed in central Ohio. Dynamic interaction was achieved by running the two models individually with an intermediate interchange of information between the surface and subsurface compartments. The main objective of this research was therefore to take advantage of TOPNET and MODFLOW, which are effective tools for detailed surface water modeling and groundwater modeling, respectively, and integrate them into a single model that can adequately simulate watershed hydrology. Models that simulate surface hydrology usually oversimplify the impact of groundwater flow processes, while groundwater models often simplify surface water flow processes. In order to overcome this simplification, there is a need for methods that can effectively simulate water flow through the unsaturated and saturated zones in large-scale hydrological models.

2 Overview of TOPNET model

MODFLOW is a standard groundwater simulation model; a detailed description can be found in McDonald and Harbaugh (1988) and will therefore not be given in this paper. TOPNET, on the other hand, is a relatively new concept in simulating rainfall-runoff processes.

A brief outline of the main components is given here. The hydrological model TOPNET (Ibbitt et al. 2001; Bandaraogoda et al. 2004) is a distributed rainfall-runoff routing model based on the TOPMODEL concepts (Beven and Kirkby 1979; Beven et al. 1995) and kinematic wave routing in a river network. The major modifications to TOPMODEL are the addition of a potential evapotranspiration (ET) component, a canopy storage component, and a soil zone component that provides the capacity of infiltration excess runoff generation using the Green-Ampt equation. There is also an additional kinematic wave channel-routing algorithm (Goring 1994). TOPNET uses a digital elevation model (DEM)-based system to delineate river channels. Most of the model components are geo-referenced and processed in ArcGIS, which is used to pre-process distributed data of watersheds and sub-watersheds and assign points of interest, such as stream gauging stations. TOPNET keeps daily accounts of the water balance components of a catchment. The water balance is monitored by observation of root zone water, which can be lost to groundwater and to ET if the groundwater level is close to the ground surface. Water in each sub-watershed flows into a river network and subsequently flows out through the end point of the watershed.

In TOPNET, all precipitation becomes either surface runoff or infiltration, according to infiltration parameters of the watershed. If the water table depth is shallow, the catchment can become saturated to the surface, resulting in saturation excess runoff generation. In addition, if the soil in the root zone is dry, then more water can infiltrate. ET is calculated by first estimating a potential ET based on temperature and day length using the Priestley-Taylor equation, and then adjusting for the increase or decrease in evaporation due to vegetation and canopy cover characteristics. If the soil in the root zone is wet enough, then the actual value of ET is the potential ET, and if the soil moisture is below the field capacity, then the actual value of ET is proportionately less than the potential value. If the soil is wet with the soil moisture in excess of the field capacity, then water drains to the shallow groundwater system. Water flows from the groundwater zone into streams as baseflow. The more water there is in the groundwater system, the faster it flows into the streams. The flow in streams is routed through the river network using kinematic wave modeling.

2.1 TOPNET model inputs

The following description of the main inputs to TOPNET is adopted from Bandaragoda et al. (2004). The main model inputs are precipitation and meteorological parameters such as wind speed and minimum and maximum temperatures. Precipitation provides input to the canopy interception, and interception storageSIis obtained by the following equation:

wherePis the precipitation rate,Cris an interception adjustment factor,Eis the reference evapotranspiration rate, andis a function providing throughfall, estimated as

whereCcis the canopy capacity. The throughfall is

Reference ET demand unsatisfied by evaporation of intercepted water is

2.2 Root zone storage component

The main parameters that describe the root zone storage processes are the depth of the root zone (d), unsaturated hydraulic conductivity (K), Green-Ampt wetting front suction (ψf), soil drainage parameter (c), drainable moisture (Δθ1), available plant moisture (Δθ2), and the impervious surface fraction (fI). Infiltration excess runoff and drainage to the saturated zone are influenced by the root zone storage. All the soil parameters except the impervious surface fraction are estimated using the Clapp and Hornberger (1978) soil textural relationships. The impervious surface fraction is determined from land cover and land use. There is no infiltration in impervious areas, and therefore infiltration is zero in these areas while surface runoff is at its maximum. In the pervious areas, the state variabledefines the amount of water held in the root zone, and it is obtained from the following equation:

whereIis the infiltration rate,Esis the soil evapotranspiration rate, andRis the soil zone drainage rate or recharge to the saturated zone.Iis assumed to be less than the infiltration capacity (Ic), which is modeled using the Green-Ampt equation:

whereZfis the depth of the wetting front and is estimated by assuming that all water in the root zone occupies a saturated zone above the wetting front. It is obtained as follows:

When there is surface water, unsatisfied ET demand is given first priority and infiltration occurs only when there is excess surface water after ET demand has been met. When this excess water exceedsIc, infiltration excess surface runoff is generated. When the moisture content is greater than the field capacity, there is drainage from the soil zone. The relative drainable saturationSrdis defined as

Recharge to the saturated zone is

whereis saturated hydraulic conductivity. Soil ET is unlimited when soil moisture content is in excess of field capacity. However, between field capacity and the permanent wilting point, ET decreases linearly, arriving at zero when the wilting point is reached.

2.3 Saturated zone component

The saturated zone component is modeled using the TOPMODEL assumptions of saturated hydraulic conductivity decreasing exponentially with soil profile depth and saturated lateral flow being driven by topographic gradients (Beven and Kirkby 1979; Beven et al. 1995). Two parameters, the soil profile lateral transmissivityand the sensitivity parameterf, characterize the decrease of hydraulic conductivity with soil profile depth. Using these TOPMODEL assumptions, a state variable, the average soil moisture deficit, is obtained as follows:

whereis the average water table depth andis the spatial average of the topographic wetness index, given by the equation

whereαis the specific catchment area, and tanβis the topographic slope. The parametersandfare estimated by means of soil textural relationships at different depths. The topographic variablesαand tanβare obtained using the TauDEM terrain analysis method developed by Tarboton (1997). As in TOPMODEL, the local water table depthzis a function of the topographic wetness index:

The distribution of wetness index is represented using a histogram of wetness index classes with the proportion of area falling within each class recorded and water table depth calculated for each class. The water table depth is used for determining areas of surface saturation, and the excess surface water input becomes saturation excess surface runoff. The water table depth in each class is also used to determine the parts of the model element where the groundwater saturated zone upwells into the soil zone, which represents loss of water from the groundwater saturated zone.

3 Study methodology

3.1 Study area and climate

The Big Darby Watershed is located 40 km west of downtown Columbus, Ohio, and covers 1 440 km2. Fig. 1 shows the location of the Big Darby Watershed, the topography, the major stream network, and the location of the two main stream gauging stations in thewatershed. The watershed covers parts of Logan, Clark, Union, Champaign, Madison, Franklin, and Pickaway counties. The general terrain of the watershed varies from rolling hills at the headwaters in Logan County, to flat plains in the middle section, to floodplains near the mouth where the Big Darby Creek meets the Scioto River. The main tributaries of the Big Darby Creek are Flat Branch, Spain Creek, Buck Run, Treacle Creek, Sugar Run, Little Darby Creek, Hellbranch Run, Spring Fork, and Robinson Run.

Fig. 1Location of Big Darby Watershed, major streams, gauging stations, and topography

The Big Darby Watershed lies in the temperate climate of central Ohio. It can be divided into three sub-watersheds: the Upper Big Darby, the Lower Big Darby, and the Little Darby sub-watersheds (Fig. 1). The Midwestern Regional Climate Center (MRCC) collects historical climate data for the Big Darby Watershed at stations located in Irwin, Marysville, and Circleville. Generally, summers are hot and humid while winters are cold and cloudy. There is usually little variation in average seasonal temperatures. Relative humidity ranges from 60% mid-afternoon to 80% in the pre-dawn hours. Average wind speed ranges from 20.9 km/h to 29.0 km/h. Thunderstorms are common from April to August. Weather data obtained from 1991 to 1997 at the Ohio weather station in Irwin are as follows: the mean monthly temperature ranged from –3℃ in January to 23.4℃ in July, with the annual mean temperature being 11.2℃; the mean monthly precipitation ranged from 49.5 mm in February to 118 mm in July, with the mean annual precipitation being 969 mm.

3.2 Model coupling methodology and governing equations

The coupled model was developed in three stages: (1) study area conceptual model and coupling model design, (2) model coupling and testing, and (3) application and evaluation. The watershed model TOPNET, a networked version of TOPMODEL, and the groundwater flow model MODFLOW-96 were selected to simulate the surface water and groundwater dynamics.

The coupling process was designed to address three major aspects: coupling components, spatial discretization and coupling, and temporal discretization and coupling. Performance of the coupled model was tested by comparing model outputs (streamflow and water table depth) with measured data for the Big Darby Watershed. As outlined by Panday and Huyakorn (2004), coupling of surface and subsurface flow models can be achieved by (1) a full coupling or fully implicit approach, (2) a sequential coupling approach in which the interaction flux is applied as a boundary condition to each model, or (3) a sequential coupling approach in which the groundwater head for one system acts as a general head boundary for the other system. This model coupling approach is based on the potential coupling interface tool developed by Bulatewicz (2006). It is a sequentially coupling approach in which the output from one model is used as the input to the other models while they run sequentially. An advantage of this sequential coupling approach is that many sub-time steps can be used for the surface model (due to the rapid surface water wave propagation speeds) before solving for the longer time steps in the subsurface flow model (Fairbanks et al. 2001).

In the coupling of TOPNET with MODFLOW, instead of measuring water table depth based on the TOPNET wetness index, the water table depth for each groundwater model cell is passed from MODFLOW to TOPNET. As a result of this modification, water table depth calculations are done for each model node, instead of the wetness index class. Recharge is an output of TOPNET and is lumped over a model element (sub-watershed).

The interactions between MODFLOW and TOPNET proceed as follows:

(1) MODFLOW provides the baseflow and water table depth at each node to TOPNET;

(2) TOPNET uses the water table depth, root zone depth, precipitation, evapotranspiration, and snowmelt to determine streamflow;

(3) TOPNET determines the net recharge to the saturated zone and passes it to MODFLOW;

(4) MODFLOW uses the recharge to calculate the water table depth.

In the coupled model, the vertical water flux from the saturated zone is calculated by TOPNET at every time step and forms the groundwater recharge to all active cells of MODFLOW. This method uses the mass conservation approach in which the leakage flux (recharge) from the surface water model is applied to the groundwater model. For each stress period in MODFLOW, hydraulic stresses are assumed to be constant.

In TOPNET, the amount of water held in the soil zone for each model element at time stepis calculated with Eq. (5). The recharge of each time step is summed up according to the following equation (Langevin et al. 2005) to obtain an average recharge estimate,Rj, for thejth stress period in MODFLOW:

whereis the recharge to the saturated zone for theith time step in thejth stress period of MODFLOW, andnis the number of time steps used in TOPNET in a stress period.

3.3 Temporal and spatial discretization

(1) Temporal discretization:

The two models use different temporal discretization schemes. Generally, surface water models use smaller time steps in the order of hours or days, while groundwater models use larger time steps. The smallest time step for groundwater modeling would be days, but use of monthly time steps is common because of the low velocities of water movement in the subsurface compared to water movement in streams. Therefore, in coupled groundwatersurface water models it is important to synchronize the time steps in order to obtain reasonable results. In this study, a monthly time step was used for the groundwater model and a daily time step for the surface water model.

(2) Spatial discretization:

An important component of the coupling of the TOPNET and MODFLOW models is the spatial linkage of the sub-watersheds used by TOPNET with the finite-difference cells used by MODFLOW. Two spatial conversions must be performed. Recharge calculated by TOPNET in a sub-watershed must be distributed over the corresponding MODFLOW cells, and water table depths for the MODFLOW cells must be combined to produce a water table elevation for each sub-watershed. GIS technology was used to join TOPNET sub-watersheds to MODFLOW grid cells by areally averaging the grid cells that fell within a particular sub-watershed. ArcMap was used to determine which MODFLOW cells fell within each sub-watershed and to areally average the water table depth in the sub-watershed for each time step.

3.4 Coupled model

Initially, the TOPNET and MODFLOW models are set up and run individually. The surface and subsurface boundaries of the coupled model are defined to be identical to those of the individual models. During the coupling development, the surface and subsurface compartments of the Big Darby Watershed are connected to each other through the three TOPNET sub-watersheds.

During the operation of the coupled model, the major rivers contribute water to the aquifers, as modeled by the River Package in MODFLOW, and the three sub-watersheds provide recharge to the groundwater, as regulated by the infiltration function of TOPNET. Both surface water and groundwater models are linked through groundwater recharge and river-aquifer interaction, using an interface program to exchange input and output parameters. The coupled parameters are transferred back and forth with the time-series output of the TOPNET-modeled percolation and sent as an input to the groundwater model. Then, the groundwater model is run again with these new input parameters and the river-groundwater interaction terms to recalculate the streamflow computed earlier with the surface water model.

The final output of the coupled model represents the results of the individual hydrological components of the original TOPNET and MODFLOW model as well as the dynamic interaction in their interface zones.

3.5 Calibration and verification

A trial and error calibration approach was used. Calibration targets were streamflow measured at the Darbyville gauging station in the Big Darby sub-watershed and at the Jefferson gauging station in the Little Darby sub-watershed, and water table depths measured with piezometers within the watershed. Daily streamflow at the two gauging stations for the period from October 1992 to September 1996 was used for model calibration. The subsequent verification of the surface water model was then performed with data from October 1996 to September 1999. Calibration and verification of MODFLOW was done for a steady state as well as a transient state. Model calibration was accomplished by varying the model input parameters within plausible ranges to produce the best fit between simulated and observed water table depths in the watershed. Water level measurements for nine wells that were used for estimation of potentiometric surfaces of the watershed aquifer were considered for calibration of the steady state model. Fig. 2 shows the location of the wells. Data from nine wells was used because of the unavailability of data for the period under consideration for the other three wells. The hydraulic head surface from the steady state simulation provided initial conditions for the transient simulation.

Fig. 2Location of monitoring wells in study area

3.6 Evaluation criteria

There are several criteria for model calibration that have been proposed and discussed (Green and Stephenson 1986; Martinec and Rango 1989; Loague and Green 1991; Refsgaard 1997; Weglarczyk 1998; Legates and McCabe 1999). A judicious combination of several techniques should be employed for a thorough model assessment. We used the Nash Sutcliff Indexgoodness of fitand percent biasto evaluate the utility of the coupled model. The performance of the model in simulating water table depths was also evaluated using the mean absolute errorstatistic.PBis a measure of the average tendency of simulated flows to be smaller or larger than the measured or observed values. Therefore, an optimum value ofPBis zero. A positivePBrepresents model underestimation, while a negative value indicates model over-prediction (Gupta et al. 1999). For calibration and validation,PBvalueshave been considered very good when they are less than 10%, good in the range 10% to 15%, and fair between 15% to 25% (Donigian et al. 1983). The same criteria were therefore adopted in this study.PBvalues greater than 25% were considered unsatisfactory. Servat and Dezetter (1991) foundENSto be the best objective function for reflecting the overall fit of a hydrograph. In this study,values greater than 0.75 were considered good.can be used to compare the model prediction and observation. The deviation offrom unity is a useful indicator of model data agreement.

4 Results and discussion

4.1 Streamflow

In order to evaluate the performance of the coupled model, daily streamflows were simulated for the period from 1992 to 1999. Fig. 3(a) shows the time series comparison of streamflow predicted using the coupled model and the measured values at the Darbyville gauging station for the period from October 1992 to September 1996. The satisfactory match can be observed in the figure, except for over-prediction which is most notable in summer of 1995. Although it is apparent that further model calibration may improve the results, over-prediction of streamflows could be attributed to a non-representative land use parameter file used in this study, due to unavailability of accurate land use maps. Errors in land use lead to inaccuracies in quantified surface runoff and ET. Simplified assumptions on aquifer hydrogeology in the groundwater model may also have resulted in over estimation of baseflow contribution to streamflow, leading to over-prediction of streamflow. However, in general, the model adequately matches the time series trend in streamflow. Close agreement between the observed streamflow and streamflow simulated with the coupled model was also achieved during the verification period from October 1996 to September 1999, shown in Fig. 3(b).

Fig. 3Measured and simulated streamflows at Darbyville gauging station

The figures show that there is underestimation of high streamflows. A possible reason could be that the effect of snow is not simulated in this model. Some other reasons for underestimationof streamflow during these high rainfall periods are the effects of localized storm events which cannot be captured by the weather stations in the watershed and probable contribution of flow from areas considered non-contributing by the model. However, the coupled model is capable of simulating the consistent overall trend for both the calibration and verification periods. A summary of statistical evaluation of streamflow for the coupled model and TOPNET is shown in Table 1.

Table 1Statistical evaluation of streamflow for calibration and validation periods

The, andENSvalues show that the coupled model is able to satisfactorily simulate streamflow in the watershed. The differences between the simulations of the coupled model and TOPNET may be due to the effect of baseflow. Quantification of baseflow in MODFLOW, which was better than quantification of baseflow in TOPNET, likely resulted in better streamflow simulation by the coupled model. Spatial and temporal changes in baseflow contributions to streamflow are mainly determined by the river-aquifer flow exchange rate in the coupled model. This exchange is affected by aquifer properties such as hydraulic conductivity, storativity, initial water table depth, and aquifer depth. However, the effects of these parameters are not captured in the TOPNET model and this affects the streamflow simulated by TOPNET. In comparison to TOPNET, the coupled model performed relatively well. For the calibration period, averageandvalues of 0.83 for the two gauging stations were obtained for the coupled model, while TOPNET yielded averageandvalues of 0.83 and 0.82 for the two gauging stations, respectively. During the validation period, the coupled model yielded averageandvalues of 0.85 and 0.84 for the two gauging stations, respectively, while TOPNET yielded 0.85 and 0.86.also showed marginal differences between the coupled model and TOPNET. An average bias of 12.9% was obtained for TOPNET during the calibration period, while the coupled model yielded a bias of 11.2% during the same period. For the validation period the coupled model yielded an average bias of 8.8%, while TOPNET yielded an average bias of 9.3%. The modest improvement in streamflow simulation using the coupled model is most likely the effect of improved baseflow simulation using MODFLOW. This modest improvement may be further enhanced if such a model is applied in a watershed where surface water-groundwater interactions are more pronounced.

4.2 Water table depth

Fig. 4 shows the contours for the mean values of simulated and measured water table depths in the wet season in 1998.

Fig. 4Contour maps of measured and simulated water table depths in wet season of 1998 in Big Darby Watershed

The spatial pattern of the contours shows how well the coupled model is able to simulate water table depths. Figs. 5(a) through (c) are representative figures showing measured and simulated water table depths in three piezometers.

Fig. 5Measured and simulated water table depths

The coupled model is able to effectively capture the general trend of groundwater dynamics. Anomalies may be due to the effect of the recharge, which is averaged for the longer time step used in MODFLOW. However, in general, the coupled model is capable of describing the average water table depth and the amplitude of its fluctuations in the Big Darby Watershed. Table 2 shows the mean absolute errorandvalues for simulated water table depths at nine monitoring wells in the watershed. Generally, the simulated depths were shallower than observed values.

Table 2Mean absolute errors and goodness of fit for simulated water table depths at nine observation wells

The influence of the tile drains in the watershed, which is not accounted for either in the coupled model or in MODFLOW, is the likely cause of this trend. The watershed is predominantly used for agriculture and there is an extensive network of irrigation tile drains. The drains were not included in the model due to lack of data on drain configuration. Table 2 shows that the mean absolute error of the simulated water table depthvaries from 0.5 m to 2.3 m for the coupled model and from 0.94 m to 2.07 m for MODFLOW. The overallfor the coupled model was 1.38 m, compared to 1.49 m for MODFLOW, during the calibration period, while the overallfor the coupled model was 1.37 m, compared to 1.50 m for MODFLOW during the validation period. From the values of goodness of fit in Table 2, we can obtain that the coupled model was superior to MODFLOW in water table depth simulation. The higher errors for MODFLOW are most likely due to the influence of groundwater recharge, which is considered to vary spatially and temporarily in the coupled model. However, water table depths simulated by the coupled model reproduce the annual variations in water table depths more accurately than MODFLOW. Although there is no consideration of spatial variation in recharge rates within a sub-watershed, the use of different recharge rates in each sub-watershed likely caused the obvious improvement in the coupled model-simulated water table depths compared to MODFLOW alone, which utilizes a single recharge rate for the entire watershed.

5 Conclusions

A model was developed by coupling TOPNET and MODFLOW using data exchange atspecified locations, and the model was applied to the Big Darby Watershed. This study evaluated the potential of using model coupling interfaces as an alternative to expensive new integrated model code writing. In this approach, data are exchanged between MODFLOW cells and TOPNET sub-watersheds (model elements). Main parameters in this approach are recharge to groundwater and water table depth. This approach replaces the wetness index-based water table depths used by TOPNET with the water table depths calculated by MODFLOW, while the average groundwater recharge per stress period used in MODFLOW is replaced with time-varying recharge simulated by TOPNET. The model was calibrated and validated with observed streamflow at two gauging stations in the watershed and measured water table depths for the period from 1992 to 1999. The coupled model was able to consistently predict annual streamflow variations at the two gauging stations, and resulted in modest improvements in streamflow simulation compared to TOPNET. Meanwhile, in comparison to MODFLOW, the coupled model consistently predicted the water table depths at the selected groundwater monitoring wells with reasonable accuracy. This research shows that the simplified model coupling approach of coordinating data transfer between models is a promising tool and can be useful whenever groundwater-surface water interaction is of concern. This study presents a methodology that can be used to assess impacts of different stresses, such as climate change and land use, on surface water and groundwater reserves. The methodology combines the advantages of a spatially distributed surface water model and a widely proven modular finite difference groundwater model. The modeling approach likely better represents the interdependency between recharge processes of surface and subsurface systems. As outlined by Goderniaux et al. (2009), using such integrated models also enables better identification of the origin of model inaccuracies in the interpretation of the results of simulations. Application of this coupled model in areas with high levels of groundwater-surface water interaction, such as wetlands and floodplains, will most likely result in much improved results of both streamflow and water table depth simulation. As stated by Markstrom et al. (2008), these coupled models should not be evaluated solely on the basis of their ability to predict streamflow at a basin outlet, but also, using different measures, on their ability to reproduce changes in surface and subsurface flows and storage in the modeled areas. Such analyses and evaluations are critical to today's scientific inquiry of and debate on sustainability of water resources.

Acknowledgements

This research was undertaken with support from the Utah Water Research Laboratory and the Department of Biological and Irrigation Engineering at Utah State University. The authors are grateful to Dr. David Tarboton for the TOPNET model code and guidance on how to use it. Many thanks to Dr. Christina Bandaragoda and Ndihui Gathuma for their help with the model setup, generation of ArcGIS soil and land use parameter files, and the use of the TOPSETUP code. The authors also thank Thomas Bulatewicz for providing expertise in model coupling.

References

Bandaragoda, C., Tarboton, D. G., and Woods, R. 2004. Application of TOPNET in the distributed model intercomparison project.Journal of Hydrology, 298, 178-201. [doi:10.1016/j.jhydrol.2004.03.038]

Beven, K. J., and Kirkby, M. J. 1979. A physically based, variable contributing area model of basin hydrology.Hydrological Sciences Journal, 24(1), 43-69. [doi:10.1080/02626667909491834]

Beven, K. J., Quinn, P., Romanowicz, R., Freer, J., Fisher, J., and Lamb, R. 1995.TOPMODEL and GRIDATB, A Users Guide to the Distribution Versions for DOS. Technical Report 110 (2nd edition). Lancaster: Centre for Research on Environmental Systems and Statistics, Lancaster University.

Bulatewicz, T. F., Jr. 2006.Support for Model Coupling: An Interface-Based Approach. Ph. D. Dissertation. Eugene: University of Oregon.

Clapp, R. B., and Hornberger, G. M. 1978. Empirical equations for some soil hydraulic properties.Water Resources Research, 14(4), 601-604.[doi:10.1029/WR014i004p00601]

Donigian, A. S., Jr., Imhoff, J. C., and Bicknell, B. R. 1983. Predicting water quality resulting from agricultural non point source pollution via simulation-HSPF. Schaller, F. W., and Baily, G. W., eds.,Agricultural Management and Water Quality, 200-249.Ames: Iowa State University Press.

Ellingson, C., and Schwartzman, P. 2004. Integration of a detailed groundwater model into a regional HSPF model.Newsletter. Golden: International Ground Water Modeling Center. http://igwmc.mines.edu/ news/spring04news.pdf [Retrieved September 2008].

Fairbanks, J., Panday, S., and Huyakorn, P. S. 2001. Comparisons of linked and fully coupled approaches to simulating conjunctive surface/subsurface flow and their interactions. Seo, B., Poeter, E., and Zheng, C., eds.,MODFLOW 2001 and Other Modeling Odysseys,Conference Proceedings, 356-361. Golden.

Goderniaux, P., Brouyère, S., Fowler, H. J., Blenkinsop, S., Therrien, R., Orban, P., and Dassargues, A. 2009. Large scale surface-subsurface hydrological model to assess climate change impacts on groundwater reserves.Journal of Hydrology, 373(1-2), 122-138. [doi:10.1016/j.jhydrol.2009.04.017]

Goring, D. G. 1994. Kinematic shocks and monoclinal waves in the Waimakariri, a steep, braided, gravel-bed river.Proceedings of the International Symposium on Waves, 336-345.Vancouver: University of British Columbia.

Green, I. R. A., and Stephenson, D. 1986. Criteria for comparison of single event models.Hydrological Sciences Journal, 31(3), 395-409. [doi:10.1080/02626668609491056]

Gupta, H. V., Sorooshian, S., and Yapo, P. O. 1999. Status of automatic calibration for hydrological models: Comparisons with multilevel expert calibration.Journal of Hydrologic Engineering, 4(2), 135-143. [doi:10.1061/(ASCE)1084-0699(1999)4:2(135)]

Ibbitt, R. P., Henderson, R. D., Copeland, J., and Wratt, D. S. 2001. Simulating mountain runoff with meso-scale weather model rainfall estimates: A New Zealand experience. Journal of Hydrology, 239(1-4), 19-32. [doi:10.1016/S0022-1694(00)00351-6]

Jones, J. P., Sudicky, E. A., Brookfield, A. E., and Park, Y. J. 2006. An assessment of the tracer-based approach to quantifying groundwater contributions to stream flow.Water Resources Research, 42(2), W02407. [doi:10.1029/2005WR004130]

Langevin, C., Swain, E., and Wolfert, M. 2005. Simulation of integrated surface-water/ground-water flow and salinity for a coastal wetland and adjacent estuary.Journal of Hydrology, 314(1-4), 212-234. [doi:10.1016/j.jhydrol.2005.04.015]

Legates, D. R., and McCabe, G. J., Jr. 1991. Evaluating use of “goodness-of-fit” measures in hydrological and hydro climatic model validation.Water Resources Research, 35(1), 233-241. [doi:10.1029/1998WR 900018]

Loague, K., and Green, R. E. 1991. Statistical and graphical methods for evaluating solute transport models: Overview and application.Journal of Contaminant Hydrology, 7(1-2), 51-73. [doi:10.1016/0169-7722(91)90038-3]

Markstrom, S. L., Niswonger, R. G., Regan, R. S., Prudic, D. E., and Barlow, P. M. 2008. GSFLOW-coupled ground-water and surface-water flow model based on the integration of the precipitation-runoff modelingsystem (PRMS) and the modular ground-water flow model (MODFLOW-2005).Ground-water/Surfacewater Book 6, Modeling Techniques. Reston: U. S. Department of the Interior, U. S. Geological Survey.

Martinec, J., and Rango, A. 1989. Merits of statistical criteria for the performance of hydrological models.Journal of the American Water Resources Association, 25(2), 421-432. [doi:10.1111/j.1752-1688. 1989.tb03079.x]

McDonald, M. G., and Harbaugh, A. W. 1988. A modular three-dimensional finite difference ground water flow model. Techniques of Water-Resources Investigations, Book 6. Reston: U. S. Department of the Interior, U. S. Geological Survey. http://pubs.usgs.gov/twri/twri6a1 [Retrieved September 2008].

Nemeth, M. S., and Solo-Gabriele, H. M. 2003. Evaluation of the use of reach transmissivity to quantify exchange between groundwater and surface water. Journal of Hydrology, 274(1-4), 145-159. [doi:10.1016/S0022-1694(02)00419-5]

Panday, S., and Huyakorn, P. S. 2004. A fully coupled physically-based spatially-distributed model for evaluating surface/subsurface flow.Advances in Water Resources, 27(4), 361-382. [doi:10.1016/j. advwatres.2004.02.016]

Pinder, G. F., and Sauer, S. P. 1971. Numerical simulation of flood wave modification due to bank storage effects.Water Resources Research, 7(1), 63-70. [doi:10.1029/WR007i001p00063]

Refsgaard, J. C. 1997. Parameterization, calibration and validation of distributed hydrological models.Journal of Hydrology, 198 (1-4), 69-97. [doi:10.1016/S0022-1694(96)03329-X]

Servat, E., and Dezetter, A. 1991. Selection of calibration of objective functions in the context of rainfall-runoff modeling in a Sudanese savannah area.Hydrological Science Journal, 36(4), 307-330. [doi:10.1080/02626669109492517]

Smits, F. C., and Hemker, C. 2004. Modeling the interaction of surface-water and groundwater flow by linking Duflow to Microflow.FEM-MODFLOW International Conference on Finite Element Models, Modflow and More. http://www.microfem.com/download/surface-grw.pdf [Retrieved 2008].

Swain, E. D., and Wexler, E. J. 1993.A Coupled Surface-water and Ground-Water Flow Model for Simulation of Stream-Aquifer Interaction, U. S. Geological Survey Open-file Report, 92-138.

Tarboton, D. G. 1997. A new method for the determination of flow directions and contributing areas in grid digital elevation models.Water Resources Research, 33(2), 309-319. [doi:10.1029/96WR03137]

Weglarczyk, S. 1998. The interdependence and applicability of some statistical quality measures of hydrological models.Journal of Hydrology, 206 (1-2), 98-103. [doi:10.1016/S0022-1694(98)00094-8]

Werner, A. D., Gallagher, M. R., and Weeks, S. W. 2006. Regional-scale, fully coupled modeling of stream-aquifer interaction in a tropical catchment.Journal of Hydrology, 328(3-4), 497-510. [doi:10.1016/j.jhydrol.2005.12.034]

*Corresponding author (e-mail:acguzha@ufl.edu)

Dec. 10, 2009; accepted Jun. 1, 2010