Simulating streamflow and water table depth with a coupled hydrological model

DEWater Science and Engineering, 2010, 3(3): 241-256doi:10.3882/.issn. 1674-2370.2010.03.001http://www.waterjournal.cne-mail: wse2008@vip.163.comSimulating streamflow and water table depth with acoupled hydrological modelAlphonce Chenjerayi GUZHA*l, Thomas Byron HARDY21. Department of Agriculural and Biological Engineering, Southwest Florida Research and Education Center,Universiy of Florida, Immokalee, FL 34142, USA2. River Systems Institute, Texas State Universiy, San Marcos, Texas 78666, USAAbstract: A coupled model integrating MODFLOW and TOPNET with the models interactingthrough the exchange of recharge and baseflow and river-aquifer interactions was developed andapplied to the Big Darby Watershed in Ohio, USA. Calibration and validation results show thatthere is generally good agreement between measured streamflow and simulated results from thecoupled model. At two gauging stations, average goodness of ft (R2 ), percent bias (P), andNash Sutcliffe eficiency ( Ems ) values of 0.83, 11.15%, and 0.83, respectively, were obtained forsimulation 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 average R2 values of 0.77and 0.76 for calibration and validation, respectively. The good match between measured andsimulated streamflows and water table depths demonstrates that the model is capable of adequatelysimulating strearmflows and water table depths in the watershed and also capturing the influence ofspatial and temporal variation in recharge.Key words: hydrological modeling; model coupling; streamflow; groundwater; TOPNET model;MODFLOW model; Big Darby Watershed1 IntroductionIn recent years, integrated surface-subsurface modeling tools have evolved rapidly andare now being applied in various watershed hydrology studies in different parts of the world.Studies include those by Markstrom et al. (2008), Panday and Huyakom (2004), Jones et al.(2006), and Wermer et al. (2006). The application of coupled models provides evidence of thecapacity of these models to produce realistic catchment behavior. However, as stated byNemeth and Solo-Gabriele (2003), linking groundwater and surface water models isfrequently problematic because the models use different sets of goverming equations.Additionally, the time scale is usually longer for groundwater modeling than for surfacewater modeling.Considerable effort has been expended to characterize the physical, chemical, and biologicalprocesses affecting groundwater and surface water resour中国煤化工pecause it has:TYHCNMHG*Corresponding author (e mail: acguzha@uf1.edu)Received Dec. 10 2009; accepted Jun. 1, 2010become 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 evaluatedmeaningfully without considering the dynamics in subsurface flow systems. As the developmentof fully integrated model concepts for this purpose is still in its early stages, one means ofintegration is the coupling of existing models. This leads to other problems, however, becausemost models have been designed to simulate specific aspects of the water cycle. Coupling of twoor more models can also result in inconsistencies because the individual models may describe thesame processes in different ways. In various studies, the coupling of surface and subsurface flowmodels has begun with the setup of relationships between river stages and groundwater storage(Pinder and Sauer 1971). Early attempts at coupling hydrological models includedMODBRANCH (Swain and Wexler 1993), which couples the groundwater flow modelMODFLOW 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 groundwatermodel 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 modelingsystem (PRMS). This approach of coupling already existing models will likely reduce costs indevelopment of new integrated model codes. In most cases where model integration has beenattempted, it is very costly to build a single predictive model that adequately represents allhydrological 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 betweengroundwater and surface water to quantify the effects of human activity and natural phenomenaon watershed hydrological responses. The research was carried out using the TOPNET andMODFLOW models with application to the Big Darby Watershed in central Ohio. Dynamicinteraction was achieved by running the two models individually with an intermediateinterchange of information between the surface and subsurface compartments. The mainobjective of this research was therefore to take advantage of TOPNET and MODFLOW, whichare 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 flowprocesses, while groundwater models often simplify surface water flow processes. In order toovercome this simplification, there is a need for methods that can effectively simulate waterflow through the unsaturated and saturated zones in large -scale hydrological models.2 Overview of TOPNET modelMODFLOW is a standard groundwater simulation model; a detailed description can befound in McDonald and Harbaugh (1988) and will therefore not he oiven in this paper.TOPNET, on the other hand, is a relatively new conce中国煤化工nof pocesse..sCNMHG242 Alphonce Chenjerayi GUZHA et al. Water Science and Engineering, Sep. 2010, Vol. 3, No. 3, 241-256A brief outline of the main components is given here. The hydrological model TOPNET (bbittet al. 2001; Bandaraogoda et al. 2004) is a distributed rainfall-runoff routing model based onthe TOPMODEL concepts (Beven and Kirkby 1979; Beven et al. 1995) and kinematic waverouting in a river network. The major modifications to TOPMODEL are the addition of apotential evapotranspiration (ET) component, a canopy storage component, and a soil zonecomponent that provides the capacity of infiltration excess runoff generation using theGreen-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 riverchannels. Most of the model components are geo-referenced and processed in ArcGIS, which isused to pre-process distributed data of watersheds and sub-watersheds and assign points ofinterest, such as stream gauging stations. TOPNET keeps daily accounts of the water balancecomponents 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 groundsurface. Water in each sub-watershed flows into a river network and subsequently flows outthrough the end point of the watershed.In TOPNET, all precipitation becomes either surface runoff or infiltration, according toinfiltration parameters of the watershed. If the water table depth is shallow, the catchment canbecome saturated to the surface, resulting in saturation excess runoff generation. In addition, ifthe soil in the root zone is dry, then more water can infiltrate. ET is calculated by firstestimating a potential ET based on temperature and day length using the Priestley-Taylorequation, and then adjusting for the increase or decrease in evaporation due to vegetation andcanopy cover characteristics. If the soil in the root zone is wet enough, then the actual value ofET is the potential ET, and if the soil moisture is below the field capacity, then the actual valueof ET is proportionately less than the potential value. If the soil is wet with the soil moisture inexcess of the field capacity, then water drains to the shallow groundwater system. Water flowsfrom the groundwater zone into streams as baseflow. The more water there is in thegroundwater system, the faster it flows into the streams. The flow in streams is routed throughthe river network using kinematic wave modeling.2.1 TOPNET model inputsThe following description of the main inputs to TOPNET is adopted from Bandaragoda etal. (2004). The main model inputs are precipitation and meteorological parameters such aswind speed and minimum and maximum temperatures. Precipitation provides input to thecanopy interception, and interception storage S is obtained by the following equation:Si=P[1-f(S)]-. EC,f(S,)(1)dtwhere P is the precipitation rate, C, is an interception adjustment factor, E is the referenceevapotranspiration rate, and f(S ) is a function provid中国煤化工STYHCNMHGAlphonce Chenjerayi GUZHA et al. Water Science and Enginering, Sep. 2010, Vol.3, No. 3, 241-256 243(S)=到(-号)(2)where C。is the canopy capacity, The throughfall isT= Pf(S.)(3)Reference ET demand unsatisfied by evaporation of intercepted water isE,=E[1-f(S)](4)2.2 Root zone storage componentThe main parameters that describe the root zone storage processes are the depth of theroot zone (d ), unsaturated hydraulic conductivity ( K ), Green-Ampt wetting front suction(Y; ), soil drainage parameter (c), drainable moisture ( 0日), available plant moisture (OO, ),and the impervious surface fraction ( f ). Infiltration excess runoff and drainage to thesaturated zone are influenced by the root zone storage. All the soil parameters except theimpervious surface ftaction are estimated using the Clapp and Hormberger (1978) soil texturalrelationships. 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 areaswhile surface runoff is at its maximum. In the pervious areas, the state variable S, definesthe amount of water held in the root zone, and it is obtained from the following equation:dS,L=I-E,-R(5)Jtwhere I is the infiltration rate, E, is the soil evapotranspiration rate, and R is the soilzone drainage rate or recharge to the saturated zone. I is assumed to be less than theinfiltration capacity (Ic ), which is modeled using the Green-Ampt equation:I.=KEr+4r(6)Zwhere Z, is the depth of the wetting front and is estimated by assuming that all water in theroot zone occupies a saturated zone above the wetting front. It is obtained as follows:Zq=-(7)△O + SO2When there is surface water, unsatified ET demand is given first priority and infiltrationoccurs only when there is excess surface water after ET demand has been met. When thisexcess water exceeds I,, infiltration excess surface runoff is generated. When the moisturecontent is greater than the field capacity, there is drainage from the soil zone. The relativedrainable saturation S。is defined asS=max(0,S, -dS02)(8)dSQRecharge to the saturated zone is中国煤化工R= K'Sq .(9)MYHCNMHG244 Alphonce Chenjerayi GUZHA et al. Water Science and Engineering, Sep. 2010, Vol. 3, No. 3, 241-256where K' is saturated hydraulic conductivity. Soil ET is unlimited when soil moisturecontent is in excess of field capacity. However, between field capacity and the permanentwilting point, ET decreases linearly, arriving at zero when the wilting point is reached.2.3 Saturated zone componentThe saturated zone component is modeled using the TOPMODEL assumptions ofsaturated hydraulic conductivity decreasing exponentially with soil profile depth and saturatedlateral flow being driven by topographic gradients (Beven and Kirkby 1979; Beven et al.1995). Two parameters, the soil profile lateral transmissivity T and the sensitivity parameterf，characterize the decrease of hydraulic conductivity with soil profile depth. Using theseTOPMODEL assumptions, a state variable, the average soil moisture deficit z△0, is obtainedas follows:d(O0)2=-R+Tge~^e~#(10)dtwhere z is the average water table depth and h is the spatial average of the topographicwetness index, given by the equationλ= ln(11)( tanB)where a is the specific catchment area, and tan B is the topographic slope. The parametersT。and f are estimated by means of soil textural relationships at diferent depths. Thetopographic variables a and tan β are obtained using the TauDEM terrain analysis methoddeveloped by Tarboton (1997). As in TOPMODEL, the local water table depth z is afunction of the topographic wetness index:_ lnI tanβ}a(12)The distribution of wetness index is represented using a histogram of wetness indexclasses with the proportion of area falling within each class recorded and water table depthcalculated for each class. The water table depth is used for determining areas of surfacesaturation, and the excess surface water input becomes saturation excess surface runoff. Thewater table depth in each class is also used to determine the parts of the model element wherethe groundwater saturated zone upwells into the soil zone, which represents loss of water fromthe groundwater saturated zone.3 Study methodology3.1 Study area and climateThe Big Darby Watershed is located 40 km west of downtown Columbus, Ohio, andcovers 1440 km'. Fig. 1 shows the location of the Big中国煤化工ography, themajor stream network, and the location of the twoYHC NM H Gations in theAlphonce Chenjerayi GUZHA et al. Water Science and Engineering. Sep. 2010, Vol. 3, No. 3, 241-256 245watershed. The watershed covers parts of Logan, Clark, Union, Champaign, Madison,Franklin, and Pickaway counties. The general terrain of the watershed varies from rolling hillsat the headwaters in Logan County, to flat plains in the middle section, to floodplains near themouth. where the Big Darby Creek meets the Scioto River. The main tributaries of the BigDarby Creek are Flat Branch, Spain Creek, Buck Run, Treacle Creek, Sugar Run, Little DarbyCreek, Hellbranch Run, Spring Fork, and Robinson Run.Ohio, Big Daibke WatershedLitle Darbyet JfferconRiver nectworkBig DarbyGauging stationElevation (m)High: 450 399DarbyilleLow 215832Fig 1 Location of Big Darby Watershed, major streams, gauging stations, and topographyThe Big Darby Watershed lies in the temperate climate of central Ohio. It can be dividedinto three sub-watersheds: the Upper Big Darby, the Lower Big Darby, and the Little Darbysub-watersheds (Fig. 1). The Midwestern Regional Climate Center (MRCC) collectshistorical 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 lttle variation in average seasonal temperatures. Relative humidity rangesfrom 60% mid-afternoon to 80% in the pre-dawn hours. Average wind speed ranges ftom20.9 km/h to 29.0 km/h. Thunderstorms are common from April to August. Weather dataobtained from 1991 to 1997 at the Ohio weather station in Irwin are as follows: the meanmonthly temperature ranged from -3C in January to 23.4C in July, with the annual meantemperature being 11.2C; the mean monthly precipitation ranged from 49.5 mm in February to118 mm in July, with the mean annual precipitation being 969 mm.3.2 Model coupling methodology and governing equationsThe coupled model was developed in three stages: (1) study area conceptual model andcoupling model design, (2) model coupling and testing, and (3) application and evaluation. Thewatershed model TOPNET, a networked version of中国煤化 indwater flowmodel MODFLOW-96 were selected to simulate theFater dynamics.YHCNMHG246 Alphonce Chenjerayi GUZHA et al. Water Science and Engineering, Sep. 2010. Vol. 3, No. 3, 241-256The coupling process was designed to address three major aspects: coupling components,spatial discretization and coupling, and temporal discretization and coupling. Performance ofthe 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 fullyimplicit approach, (2) a sequential coupling approach in which the interaction flux is applied asa boundary condition to each model, or (3) a sequential coupling approach in which thegroundwater head for one system acts as a general head boundary for the other system. Thismodel coupling approach is based on the potential coupling interface tool developed byBulatewicz (2006). It is a sequentially coupling approach in which the output from one modelis used as the input to the other models while they run sequentially. An advantage of thissequential 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 timesteps in the subsurface flow model (Fairbanks et al. 2001).In the coupling of TOPNET with MODFLOW, instead of measuring water table depthbased on the TOPNET wetness index, the water table depth for each groundwater model cell ispassed from MODFLOW to TOPNET. As a result of this modification, water table depthcalculations are done for each model node, instead of the wetness index class. Recharge is anoutput ofTOPNET 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 byTOPNET at every time step and forms the groundwater recharge to all active cells ofMODFLOW. This method uses the mass conservation approach in which the leakage flux(recharge) ftom the surface water model is applied to the groundwater model. For each stressperiod 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 stepO; is calculated with Eq. (5). The recharge of each time step is summed up according to thefollowing equation (Langevin et al. 2005) to obtain an average recharge estimate, R, for thejth stress period in MODFLOW:之oR,R,=互中国煤化工(13)iMYHCNMHGAlphonce Chenjerayi GUZHA et al. Water Science and Engineering, Sep. 2010 Vol. 3, No. 3, 241-256 247where R，is the recharge to the saturated zone for the ith time step in the jth stress period ofMODFLOW, and n is 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 watermodels use smaller time steps in the order of hours or days, while groundwater models uselarger time steps. The smallest time step for groundwater modeling would be days, but use ofmonthly time steps is common because of the low velocities of water movement in thesubsurface compared to water movement in streams. Therefore, in coupled groundwater-surface water models it is important to synchronize the time steps in order to obtain reasonableresults. In this study, a monthly time step was used for the groundwater model and a daily timestep for the surface water model.(2) Spatial discretization:An important component of the coupling of the TOPNET and MODFLOW models is thespatial linkage of the sub-watersheds used by TOPNET with the finite difference cells used byMODFLOW. Two spatial conversions must be performed. Recharge calculated by TOPNET in asub-watershed must be distributed over the corresponding MODFLOW cells, and water tabledepths for the MODFLOW cells must be combined to produce a water table elevation for eachsub-watershed. GIS technology was used to join TOPNET sub-watersheds to MODFLOW gridcells by areally averaging the grid cells that fell within a particular sub-watershed. ArcMap wasused to determine which MODFLOW cells fell within each sub-watershed and to areally averagethe water table depth in the sub-watershed for each time step.3.4 Coupled modelInitially, the TOPNET and MODFLOW models are set up and run individually. Thesurface and subsurface boundaries of the coupled model are defined to be identical to those ofthe individual models. During the coupling development, the surface and subsurfacecompartments of the Big Darby Watershed are connected to each other through the threeTOPNET sub-watersheds.During the operation of the coupled model, the major rivers contribute water to theaquifers, as modeled by the River Package in MODFLOW, and the three sub-watershedsprovide recharge to the groundwater, as regulated by the infiltration function of TOPNET. Bothsurface water and groundwater models are linked through groundwater recharge andriver-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 theTOPNET-modeled percolation and sent as an input-lel. Then, the中国煤化工groundwater model is run again with these new inpuTCHCNMHGr-groundwaterinteraction terms to recalculate the streamflow compucu C.. I . u Hi u.... water model.248 Alphonce Chenjerayi GUZHAet al. Water Science and Engineering, Sep. 2010, Vol. 3, No. 3, 241-256The final output of the coupled model represents the results of the individual hydrologicalcomponents of the original TOPNET and MODFLOW model as well as the dynamicinteraction in their interface zones.3.5 Calibration and verificationA trial and error calibration approach was used. Calibration targets were streamflowmeasured at the Darbyville gauging station in the Big Darby sub-watershed and at the Jeffersongauging station in the Little Darby sub-watershed, and water table depths measured withpiezometers within the watershed. Daily streamflow at the two gauging stations for the periodfrom October 1992 to September 1996 was used for model calibration. The subsequentverification of the surface water model was then performed with data from October 1996 toSeptember 1999. Calibration and verification ofMODFLOW was done for a steady state as well astransient state. Model calibration was+accomplished by varying the modelinput04parameters within plausible ranges to produce thebest fit between simulated and observed waterable depths in the watershed. Water levelFRIOmeasurements for nine wells that were used forM2estimation of potentiometric surfaces of theM4watershed aquifer were considered for calibrationR18of the steady state model. Fig. 2 shows the locationof the wells. Data from nine wells was usedbecause of the unavailability of data for the period●Monitoring wellunder consideration for the other three wells. Thehydraulic head surface from the steady stateFig. 2 Location of monitoring wells insimulation provided initial conditions for thestudy areatransient simulation.3.6 Evaluation criteriaThere are several criteria for model calibration that have been proposed and discussed(Green and Stephenson 1986; Martinec and Rango 1989; Loague and Green 1991; Refsgaard1997; Weglarczyk 1998; Legates and McCabe 1999). A judicious combination of severaltechniques should be employed for a thorough model assessment. We used the Nash SutcliffIndex ENs, goodness of fit R2 , and percent bias P to evaluate the utility of the coupledmodel. The performance of the model in simulating water table depths was also evaluated usingthe mean absolute error ( Az ) statistic. Pg is a measure of the average tendency of simulatedflows to be smaller or larger than the measured or observed values. Therefore, an optimumvalue of P is zero. A positive Pg represents modelYH中国煤化工egative valueindicates model over-prediction. (Gupta et al.1999). FCN MH Gn, Pg valuesAphonce Chenjerayi GUZHA et al. Water Science and Engineering, Sep. 2010, Vol. 3, No. 3, 241-256 249have 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 adoptedin this study. Pg values greater than 25% were considered unsatisfactory. Servat and Dezetter(1991) found ENs to be the best objective function for reflecting the overall fit of ahydrograph. In this study, Ens values greater than 0.75 were considered good. R2 can beused to compare the model prediction and observation. The deviation of R2 from unity is auseful indicator of model data agreement.4 Results and discussion4.1 StreamflowIn order to evaluate the performance of the coupled model, daily streamflows weresimulated for the period from 1992 to 1999. Fig. 3(a) shows the time series comparison ofstreamflow predicted using the coupled model and the measured values at the Darbyvillegauging station for the period from October 1992 to September 1996. The satisfactory matchcan be observed in the figure, except for over-prediction which is most notable in summer of1995. Although it is apparent that further model calibration may improve the results,over-prediction of streamflows could be atributed to a non-representative land useparameter file used in this study, due to unavailability of accurate land use maps. Errors inland use lead to inaccuracies in quantified surface runoff and ET. Simplified assumptions onaquifer hydrogeology in the groundwater model may also have resulted in over estimation ofbaseflow contribution to streamflow, leading to over-prediction of streamflow. However, ingeneral, the model adequately matches the time series trend in streamflow. Close agreementbetween the observed streamflow and streamflow simulated with the coupled model wasalso achieved during the verification period from October 1996 to September 1999, shownin Fig. 3(b).--. Measured- Simulated600r1500s00400告100030020000 t10Time(a) From October 1992 to Sepcember 1996(b) From October 1996 to September 1999Fig. 3 Measured and simulated streamflows at Darbyville gauging stationThe figures show that there is underestimation of k: reason couldbe that the effect of snow is not simulated in this model.中国煤化Ideresintion^TYHCNMHG250 Alphonce Chenjerayi GUZHA et al. Water Science and Engineering, Sep. 2010, Vol. 3, No. 3, 241-256of streamflow during these high rainfall periods are the effects of localized storm events whichcannot be captured by the weather stations in the watershed and probable contribution of flowfrom areas considered non-contributing by the model. However, the coupled model is capable ofsimulating the consistent overall trend for both the calibration and verification periods. Asummary of statistical evaluation of streamflow for the coupled model and TOPNET is shown inTable 1.Table 1 Statstical evaluation of streamflow for calibration and validation periodsQ (m2/s)R”p (%)EsPeriodGauging station2。2-lcRR_PERSCalibrationDarbyille25.5020.219.60.85 0.8311.6 .12.20.85 0.81(Oel.1992-SepWest Jeferson12.109.610.20.78 0.82 14.010.10.81 0.841996)Validation61.8057.45s.50.890.8110.913.30.86 0.841996-Sep.West Jefferson30.2022.824.70.830.867.:.21999)Note: %，&, and & are mean daily streamflow Q obtained from observation, TOPNET, and the coupled model,respectively;肝and 陉are average goodness of fit for TOPNET and the coupled model; PT and Pf are percent bias forTOPNET and the coupled model; and ER and路are Nash Sulcliffe eficieney for TOPNET and the coupled model.The P，R2, and Es values show that the coupled model is able to satisfactorilysimulate streamflow in the watershed. The differences between the simulations of thecoupled model and TOPNET may be due to the effect of baseflow. Quantification obaseflow in MODFLOW, which was better than quantification of baseflow in TOPNET,likely resulted in better streamflow simulation by the coupled model. Spatial and temporalchanges in baseflow contributions to streamflow are mainly determined by the river-aquiferflow exchange rate in the coupled model. This exchange is affected by aquifer propertiessuch 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 thisaffects the streamflow simulated by TOPNET. In comparison to TOPNET, the coupledmodel performed relatively well. For the calibration period, average Ens and R2 values of0.83 for the two gauging stations were obtained for the coupled model, while TOPNETyielded average ENs and R' values of 0.83 and 0.82 for the two gauging stations,respectively. During the validation period, the coupled model yielded average ENs and R2values of 0.85 and 0.84 for the two gauging stations, respectively, while TOPNET yielded0.85 and 0.86. Pg also showed marginal differences between the coupled model andTOPNET. 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 validationperiod the coupled model yielded an average bias of 8.8%, while TOPNET yielded anaverage bias of 9.3%. The modest improvement in streamflow simulation using the coupledmodel is most likely the effect of improved baseflow simulation using MODFLOW. Thismodest improvement may be further enhanced if su中国煤化工a watershedwhere surface water-groundwater interactions are mor:YHCNMHGAlphonce Chenjerayi GUZHA et al. Water Science and Engineering, Sep. 2010, Vol. 3, No. 3, 241-256 2514.2 Water table depthFig. 4 shows the contours for the mean values of simulated and measured water tabledepths in the wet season in 1998.Water table deph (m)Water lable depth (m)High. 861High: 962Low: 0.30, Low 090(a) Measured watler table depth distributon(b) Simulated water table depth distributionFig. 4 Contour maps of measured and simulated water table depths in wet season of 1998 inBig Datby WatershedThe spatial patterm of the contours shows how well the coupled model is able to simulatewater table depths. Figs. 5(a) through (c) are representative figures showing measured andsimulated water table depths in three piezometers.一Measured --- Simulated喜651 Mw多3Time@) Moutoring well LO3(b) Monitoring well M2(C) Moritoring well PK9中国煤化工Fig. 5 Measured and simulated w.MYHCNMHG252 Alphonce Chenjerayi GUZHA et al. Water Science and Engineering, Sep.2010, Vol. 3, No. 3, 241-256The coupled model is able to effectively capture the general trend of groundwaterdynamics. Anomalies may be due to the effect of the recharge, which is averaged for the longertime step used in MODFLOW. However, in general, the coupled model is capable of describingthe average water table depth and the amplitude of its fluctuations in the Big Darby Watershed.Table 2 shows the mean absolute error△z and R2 values for simulated water table depths atnine monitoring wells in the watershed. Generally, the simulated depths were shallower thanobserved values.Table 2 Mean absolute errors and goodness of ft for simulated water table depthsat nine observation wells正(m)_p?StagePeriodCoupled modelMODFLOWOct. 1992- Sep. 19931.580.760.74Oet. 1993-Sep. 1994).50.940.810.75CalibrationOct.1994 Sep.1995.61.650.770.71Oct. 199-Sep. 19962.31.770.68Oct. 1996-Sep. 1997.71.280.67ValidationOct. 1997-Sep. 1998.91.160.78Oct. 1998 Sep.1999.52.070.69The influence of the tile drains in the watershed, which is not accounted for either in thecoupled model or in MODFLOW, is the likely cause of this trend. The watershed ispredominantly 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 2shows that the mean absolute error of the simulated water table depth Az varies from 0.5 m to2.3 m for the coupled model and from 0.94 m to 2.07 m for MODFLOW. The overall Qz forthe coupled model was 1.38 m, compared to 1.49 m for MODFLOW, during the calibrationperiod, while the overall △z for the coupled model was 1.37 m, compared to 1.50 m forMODFLOW during the validation period. From the values of goodness of ft in Table 2, we canobtain 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 groundwaterrecharge, which is considered to vary spatially and temporarily in the coupled model.However, water table depths simulated by the coupled model reproduce the annual variationsin water table depths more accurately than MODFLOW. Although there is no considerationof spatial variation in recharge rates within a sub-watershed, the use of different rechargerates in each sub-watershed likely caused the obvious improvement in the coupledmodel-simulated water table depths compared to MODFLOW alone, which utilizes a singlerecharge rate for the entire watershed.5 Conclusions中国煤化工A model was developed by coupling TOPNET.fYHCNMHGexchange atAlphonce Chenjerayi GUZHA et al. Water Science and Engineering, Sep 2010, Vol. 3, No. 3, 241-256 253specified locations, and the model was applied to the Big Darby Watershed. This studyevaluated the potential of using model coupling interfaces as an altermative to expensive newintegrated model code writing. In this approach, data are exchanged between MODFLOW cellsand TOPNET sub-watersheds (model elements). Main parameters in this approach are rechargeto groundwater and water table depth. This approach replaces the wetness index-based watertable depths used by TOPNET with the water table depths calculated by MODFLOW, while theaverage groundwater recharge per stress period used in MODFLOW is replaced withtime-varying recharge simulated by TOPNET. The model was calibrated and validated withobserved streamflow at two gauging stations in the watershed and measured water table depthsfor the period from 1992 to 1999. The coupled model was able to consistently predict annualstreamflow variations at the two gauging stations, and resulted in modest improvements instreamflow simulation compared to TOPNET. Meanwhile, in comparison to MODFLOW, thecoupled model consistently predicted the water table depths at the selected groundwatermonitoring wells with reasonable accuracy. This research shows that the simplified modelcoupling approach of coordinating data transfer between models is a promising tool and can beuseful whenever groundwater surface water interaction is of concern. This study presents amethodology that can be used to assess impacts of different stresses, such as climate changeand land use, on surface water and groundwater reserves. The methodology combines theadvantages of a spatially distributed surface water model and a widely proven modular finitedifference groundwater model. The modeling approach likely better represents theinterdependency between recharge processes of surface and subsurface systems. As outlined byGodermniaux et al. (2009), using such integrated models also enables better identification of theorigin of model inaccuracies in the interpretation of the results of simulations. Application ofthis coupled model in areas with high levels of groundwater-uface water interaction, such aswetlands and floodplains, will most likely result in much improved results of both streamflowand water table depth simulation. As stated by Markstrom et al. (2008), these coupled modelsshould not be evaluated solely on the basis of their ability to predict streamflow at a basinoutlet, but also, using different measures, on their ability to reproduce changes in surface andsubsurface flows and storage in the modeled areas. Such analyses and evaluations are critical totoday's scientific inquiry of and debate on sustainability of water resources.AcknowledgementsThis research was undertaken with support from the Utah Water Research Laboratoryand the Department of Biological and Irrigation Engineering at Utah State University. Theauthors are grateful to Dr. David Tarboton for the TOPNET model code and guidance onhow to use it. Many thanks to Dr. Christina Bandaragoda and Ndihui Gathuma for their helpwith the model setup, generation of ArcGIS soil and land use parameter files, and the use ofthe TOPSETUP code. The authors also thank Thoma中国煤化宝xperise inmodel coupling.TYHCNMHG254 Alphonce Chenjerayi GUZHA et al. Water Science and Engineering, Sep. 2010, Vol. 3, No. 3, 241-256ReferencesBandaragoda, C, Tarboton, D. G, and Woods, R. 2004. Application of TOPNET in the distributed modelintercomparison project. Journal of Hydrology, 298, 178 201. [doi:10.016/jjhydrol.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:l0. 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 Hormberger, G M. 1978. Empirical equations for some soil hydraulic properties. WaterResources Research, 14(4), 601- 604.[doi: 10.1029/W R014i004p00601]Donigian, A. S, J.r, Imhoff, J. C, and Bicknell, B. R. 1983. Predicting water quality resulting fromagricultural non point source pollution via simulation-HSPF. Schaller, F. W., and Baily, G W., eds,Agricultural Management and Water Qualiy, 200-249. Ames: lowa State University Press.Ellingson, C., and Schwartzman, P. 2004. Integration of a detailed groundwater model into a regional HSPFmodel. Newsletter. Golden: Intemational Ground Water Modeling Center. htp://igwmc.mines.edu/news/spring04news.pdf [Retrieved September 2008].Fairbanks, J, Panday, S., and Huyakom, P. S. 2001. Comparisons of linked and fully coupled approaches tosimulating 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, Brouyere, s.. Fowler, H. J, Blenkinsop, S, Therrien, R., Orban, P., and Dassargues, A.2009.Large scale surfacc-subsurface hydrological model to assess climate change impacts on groundwaterreserves. Journal of Hydrology, 373( 1-2), 122-138. [doi:10.1016/jjhydrol.2009.04.017]Goring, D. G 1994. Kinematic shocks and monoclinal waves in the Waimakariri, a steep, braided, gravel-bedriver. Proceedings of the International Symposium on Waves, 336-345. Vancouver: University of BritishColumbia.Green, 1. R. A., and Stephenson, D. 1986. Criteria for comparison of single event models. HydrologicalSciences 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.106 1/(ASCE)1084-0699(1999)4:2(135)]Ibbitt, R. P., Henderson, R. D., Copeland, J, and Wratt, D. S. 2001. Simulating mountain runoff withmeso-scale weather model rainfall estimates: A New Zealand experience. Joumal of Hydrology, 239(1-4),19-32. [doi:10.1016/0022-1694000351-6]Jones, J. P, Sudicky, E. A., Brookfield, A. E., and Park, Y. J 2006. An asessment of the tracer-based approachto 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 andsalinity for a coastal wetland and adjacent estuary. Journal of Hydrology, 314(1-4)， 212-234.[doi:10.101 6/j.jhydrol.2005.04.015]Legates, D. R.. and McCabe, G J., Jr.1991. Evaluating use of "goodness-of fit measures in hydrological andhydro climatic model validation. Water Resources Research, 35(1), 233-241. [doi:10.1029/1998WR900018]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中国煤化工:LOW-coupledground-water and surface-water flow model based on the inMYHCNMHG_Inoff modelingAlphonce Chenjerayi GUZHA et al. Water Science and Engineering, Sep. 2010, Vol. 3, No.3, 241-256 255system (PRMS) and the modular ground-water flow model (MODFLOW 2005). Ground-water/Surface-water 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 performnance of hydrological models.Jourmal of the American Water Resources Association, 25(2), 421-432. [oi.10111j.1752-16881989.tb03079.x]McDonald, M. G, and Harbaugh, A. W. 1988. A modular three dimensional finite dfference ground waterflow model. Techniques of Water-Resources Investigations, Book 6. Reston: U. S. Department of theInterior, U. S. Geological Survey. htp://pubsusgs.gov/twri/wri6al [Retrieved September 2008].Nemeth, M. S, and Solo-Gabriele, H. M.2003. Evaluation of the use of reach transmissivity to quantifyexchange between groundwater and surface water. Journal of Hydrology, 274(1-4), 145-159.{doi: 10.1016/S0022-1 694(02)00419-5]Panday, S., and Huyakom, P. s. 2004. A fully coupled physically-based spatially-distributed model forevaluating 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 storageeffects. Water Resources Research, 7(1), 63-70. [doi:10.1029/WR007i001p00063]Refsgaard, J. C.1997. Parameterization, calibration and validation of ditributed hydrological models. Journalof 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 ofrainfall-runoff modeling in a Sudanese savannah area. Hydrological Science Joumal, 36(4), 307-330.[do: 10.1080/02626669109492517]Smits, F. C, and Hemker, C.2004. Modeling the interaction of surface-water and groundwater flow by linkingDuflow to Microflow. FEM-MODFLOW International Conference on Finite Element Models, Modflowand More. htp://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 Simulationof 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 griddigital 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 ofhydrological models. Journal of Hydrology, 206 (1-2), 98- 103. [doi:10.1016/S0022-1694(98)00094-8]Wermner, A. D., Gallagher, M. R., and Weeks, S. W. 2006. Regional-scale, fully coupled modeling ofstream-aquifer interaction in a tropical catchment. Journal of Hydrology, 328(3-4), 497-510.[doi: 10.1016/jjhydrol.2005.12.034]中国煤化工MHCNMHG256 Alphonce Chenjerayi GUZHA et al. Water Science and Engineering, Sep. 2010, Vol. 3, No. 3, 241-256