## Abstract

To analyze the water budget under human influences in the Huaihe River plain region in China, the authors have developed a numerical modeling system that integrates water flux algorithms into a platform created by coupling a soil moisture model with the modular three-dimensional finite-difference groundwater flow model (MODFLOW). The modeling system is largely based on physical laws and employs a numerical method of the finite difference to simulate water movement and fluxes in a horizontally discretized watershed or field. The majority of model parameters carry physical significance and can be determined by field and laboratory measurements or derived from watershed characteristics contained in GIS and remote sensing data. Several other empirical parameters need to be estimated by model calibration. The numerical modeling system is calibrated in the Linhuanji catchment (2 560 km^{2}) to estimate surface runoff, groundwater recharge, and groundwater loss for evapotranspiration and stream baseflow. Model validation is conducted at a small runoff experimental field (1.36 km^{2}) in the Wuduogou Hydrological Experimental Station to test the model’s capability to simulate hydrological components and estimate water fluxes using observed stream stage and groundwater data, as well as lysimeter-measured precipitation recharge and groundwater loss. As proven by the promising results of model testing, this physically based and distributed-parameter model is a valuable contribution to the ever-advancing technology of hydrological modeling and water resources assessment.

## 1. Introduction

In the highly populated Huaihe River plain region in China (Fig. 1), human activities, especially agriculture, influence the hydrological processes, and water resources management must focus on sustainability, tracked based on an accurate understanding of water distribution and fluxes. In other words, water storage in and movement among all dynamically linked reservoirs must be estimated in order to evaluate water availability caused by human impact. This task of modeling the water transfer among the atmosphere, watershed surface, soil column, and groundwater aquifers has provided hydrologists many challenges in understanding and evaluating the dynamic interactions of these reservoirs (National Research Council 2004). Scientists and engineers in different disciplines often focus on certain parts of the hydrological system and treat other water balance components in a more simplified manner. In engineering and earth science investigations there has been a long tradition of treating soil moisture below the root zone as groundwater, which is used as a “boundary condition” in the simulation of soil moisture dynamics. In river hydraulics and hydrodynamics of open channels, a porous subsurface is seldom considered to be an active participant of in-channel processes and dynamics. In atmospheric science, soil moisture and groundwater are represented as “buckets” of limited sizes in which water movement is not coupled with streamflow dynamics. For other scientists and resource managers, groundwater is represented as an infinitely large reservoir in which the slow-flow process is unlikely to play a significant role in the hydrological cycle at the time scales of regular human activities (Duffy et al. 2006). As a result of the strong emphasis by water resources managers and decision makers recently put on integrated (i.e., large scale) and sustainable (i.e., long term) water resource management, integrated assessment and modeling techniques have increasingly gained in popularity since the mid-1990s (Fedra and Jamieson 1996; Jamieson and Fedra 1996a, b; Dunn et al. 1996; Reitsma 1996; Andreu et al. 1996; Parker et al. 2002).

Traditionally, hydrologists developed and applied conceptual hydrological models to simulate the rainfall–runoff relationship based on simplified and empirical descriptions of the runoff generation processes [see an excellent summary by Singh (1995)]. Because of the lack of data on physical basin characteristics and the limitations of computing power, conceptual and often lumped-parameter models mainly focus on the simulation of aggregated output (total streamflow) at the watershed outlet and usually cannot offer sufficient details and accurate estimation of water fluxes in a spatially heterogeneous domain. Although these watershed hydrological models have been widely and successfully used in flood forecasting and regional water resources planning (Shen 1992; Guo et al. 1997; Xu and Guo 1994), they cannot simulate the physical processes of water transfer in a spatially discretized system, and thus offer us a quantitative evaluation of water dynamics.

As a result of technological advancements, development of physically based and distributed-parameter models have dominated the field over the past two decades. Several commercially available products such as the Danish Hydraulic Institute’s European Hydrological System’s derivative MIKE SHE (Refsgaard and Storm 1995), the Swiss Federal Institute of Technology’s Water Balance Simulation Model (WaSiM-ETH) (Schulla and Jasper 2001), and HydroGeoSphere (Therrien et al. 2004) have become internationally well known and gained popularity for use in applications to solve many kinds of water resources problems. MIKE SHE and WaSiM-ETH dynamically link the unsaturated zone model (Richard’s equation) and a 2D or 3D groundwater model in a numerical method that can simulate the entire land phase of the hydrologic cycle. MIKE SHE has become a more popular tool and has a proven track record in numerous consultancy and research applications around the world. HydroGeoSphere is a fully integrated 3D model that can simulate water flow, heat flow, and advective–dispersive solute transport on a 2D land surface and in a 3D subsurface domain under variably saturated, heterogeneous geologic conditions. Full coupling of surface and subsurface flow regimes is accomplished implicitly by simultaneously solving one system of nonlinear discrete equations describing flow and transport in both flow regimes. Obviously these full-fledged and well-packaged models can offer powerful modeling capabilities for sophisticated simulation of watershed hydrological processes. However, their applicability in a large watershed is sometimes seriously limited because physical parameters of basin characteristics are simply not available. Furthermore, hydrologists may find it difficult to modify or expand these models in order to build a comprehensive tool for water resource management. Since our ultimate goal is to develop such a tool for Huaihe River Water Resources Commission, in this study we found it impossible to use these existing models.

Soil layers and groundwater aquifers are the two main regulators of water movement and fluxes. Given the background and motivation of this study, we formulated our modeling strategies and adopted an approach to integrating a soil moisture model with a groundwater model for numerical simulation of the hydrological cycle in a vertically and horizontally discretized domain at the watershed or field scale. Many researchers have used numerical methods to solve the Richard’s equation and thus construct soil hydrological models for 1D, 2D, or 3D soil moisture simulation (Lappala et al. 1987; Ross 1990; Stothoff 1995; Simunek et al. 1998; Hsieh et al. 2000; Fayer 2000). In a plain region, soil moisture movement in the vertical direction plays a dominant role and therefore a 1D model should be sufficient. Chen and Hu (2004) developed such a model to investigate groundwater influences on soil moisture and surface evaporation. This model was adopted in the present study to simulate the downward (infiltration and percolation) and upward (evapotranspiration) movement of soil moisture. Since these water fluxes interact with groundwater dynamics, a groundwater model must be integrated with the soil moisture model. The modular three-dimensional finite-difference groundwater flow model (MODFLOW), a widely used groundwater model of the U.S. Geological Survey, has a modular structure and open environment (McDonald and Harbaugh 1988). It was selected for the model integration in order to link its three packages for simulating water exchanges between aquifer and other hydrological components (i.e., recharge, evapotranspiration, and river) with the soil moisture model and streamflow routing algorithms.

The modeling system was calibrated in the Linhuanji catchment of the Huaihe River basin to estimate surface runoff, groundwater recharge, and groundwater losses for evapotranspiration and stream baseflow. It was further tested at the Wudaogou Hydrological Experimental Station for validating the model’s capability in simulation of hydrological components and estimation of water fluxes using observed stream water stage and groundwater data, as well as lysimeter-measured precipitation recharge and groundwater loss. Model testing and application have demonstrated its accuracy and usefulness for quantitative assessment of the dynamics of water movement and storage, which is controlled by both meteorological conditions and agricultural land use in the plain region of the Huaihe River basin.

## 2. Model components and algorithms

### a. Soil moisture model

Soil moisture dynamics driven by climate fluctuations play a key role in the simulation of water transfer among ground surface, unsaturated zones, and aquifer. Soil moisture variation in the model is described by Richard’s equation. Integrating Richard’s equation through four soil layers under the assumption of vertically homogeneous soil hydraulic properties within each layer yields the following equations (Chen and Dudhia 2001; Chen and Hu 2004):

and

where subscript *i* = 1, 2, 3, and 4 is soil layer index; *d _{i}* is thickness of the

*i*th soil layer;

*θ*is soil moisture content;

_{i}*P*is precipitation

_{d}*P*falling on the ground;

*I*is the irrigated water from stream and aquifer;

_{r}*R*is surface runoff;

_{s}*K*is vertical unsaturated soil hydraulic conductivity; and

_{i}*D*the soil water diffusivity. Both

*K*and

*D*are functions of soil moisture content

*θ*and are computed from

*K*(

*θ*) =

*K*(

_{s}*θ*/

*θ*)

_{s}^{2b+3}and

*D*(

*θ*) =

*K*(

*θ*)(∂Ψ/∂

*θ*), where

*ψ*is soil water tension function and Ψ(

*θ*) = Ψ

_{s}/(

*θ*/

*θ*)

_{s}^{b}in which b is a curve-fitting parameter (Cosby et al. 1984);

*E*

_{dir}is the evaporation from the top soil surface, and

*E*in Eqs. (1)–(3) is the transpiration by vegetation through roots. The groundwater effect on soil moisture is taken into account in Eq. (4) by allowing upward soil moisture exchange between the deepest model soil layer and the groundwater table [

_{Ti}*D*(∂

*θ*/∂

*z*)

_{4}] (Chen and Hu 2004). The vertical derivative of

*D*(∂

*θ*/∂

*z*)

_{i}is approximated by

*D*[(

*θ*

_{i+1}−

*θ*)/

_{i}*d*].

_{i}### b. Groundwater model

Interactions between unsaturated and saturated zones (e.g., the groundwater recharge and loss, and between stream water and groundwater including baseflow), depend on the groundwater dynamics, which is described by the governing equation:

where *h* is the hydraulic head; *K _{x}*,

*K*, and

_{y}*K*are the hydraulic conductivities along the

_{z}*x*,

*y*, and

*z*axes;

*S*is the specific storage of the aquifer (specific yield

_{s}*S*divided by the aquifer thickness); and

_{y}*W*is a volumetric flux per unit area representing sources and/or sinks of water. It includes groundwater recharge from soil moisture

_{s}*P*

_{rg}and groundwater loss through capillary flow upward

*G*

_{up}, water exchanges between stream channel and aquifer

*Q*, and groundwater pumping rate

_{g}*W*. In MODFLOW, this equation is numerically solved using a finite-difference method and the FORTRAN source codes of MODFLOW-2000 can be easily modified for model integration and other customized applications (McDonald and Harbaugh 1988; Harbaugh et al. 2000).

_{g}### c. Simulation of water fluxes

The two models can simulate the water movement and dynamics within the unsaturated and saturated zones, respectively. The key for integrating the two models is to simulate the water fluxes, which link the two storage reservoirs together. The entire hydrological cycle involves water fluxes at three major interfaces: 1) evapotranspiration and surface runoff at the interface between atmosphere and soils; 2) groundwater recharge and loss at the interface between unsaturated and saturated zones; and 3) water exchanges (baseflow or groundwater recharge) at the interface between aquifer and stream channel.

#### 1) Evapotranspiration

In the soil moisture model, total evapotranspiration (ET) is the sum of 1) direct evaporation from the top shallow soil layer, *E*_{dir}; 2) transpiration via canopy and roots, *E _{T}*; and 3) evaporation of precipitation intercepted by the canopy,

*E*.

_{c}A simple linear method is used to calculate *E*_{dir} (Mahfouf and Noilhan 1991):

where *σ _{f}* is the green vegetation fraction (cover);

*β*= (

*θ*

_{1}−

*θ*/

_{w}*θ*

_{ref}−

*θ*), in which

_{w}*θ*

_{ref}and

*θ*is field capacity and wilting point, respectively;

_{w}*E*is the potential evaporation calculated using a Penman-based energy balance approach that includes a stability-dependent aerodynamic resistance (Mahrt and Ek 1984). Here

_{p}*E*is calculated by

_{T}where *B _{c}* is a function of canopy resistance,

*W*is intercepted canopy water content estimated from the budget for intercepted canopy water,

_{c}*S*is the maximum canopy capacity, and

*n*= 0.5. Finally, the third component of ET,

*E*, can be estimated by

_{c}The budget for intercepted canopy water is

where *P* is total precipitation. If *W _{c}* exceeds

*S*, the excess precipitation or drip,

*D*, reaches the ground [note that

_{p}*P*= (

_{d}*1−σ*)

_{f}*P*+

*D*in Eq. (1)].

_{p}#### 2) Surface runoff

In the semihumid region of China, infiltration excess overland flow and saturation overland flow can be generated from precipitation. The former surface runoff, *R _{s}*, is defined as the excess of precipitation, which does not infiltrate into the soil (

*R*−

_{s}= P_{d}*I*

_{max}). The maximum infiltration,

*I*

_{max}, is formulated as

where *K*_{1} is the upper-layer soil hydraulic conductivity and *I _{f}* is the infiltration capacity related to precipitation intensity, soil moisture deficit, and rainfall duration (Chen and Dudhia 2001).

In the wet season, the upper layer of soil may become saturated during a rainfall event, resulting in overland flow (*R _{s}* = max{

*P*−

_{d}*D*

_{x1}, 0}, where

*D*

_{x1}is the upper-layer soil moisture deficit per unit area in a simulation time step).

#### 3) Groundwater recharge and loss

Water flux that crosses the interface between saturated and unsaturated zones is either groundwater recharge from soil moisture driven by gravity or groundwater loss in soil layers driven by capillary force. The water flux *W _{e}* can be estimated by the following equation:

where *D*((∂*θ*/∂*z*))_{4} = *D*(*θ*_{4} − *θ _{s}*/

*Z*), and

_{g}*Z*is the distance between the groundwater table and the midpoint of the soil layer located immediately above the groundwater table. A positive

_{g}*W*represents groundwater recharge

_{e}*P*

_{rg}and a negative

*W*represents groundwater loss

_{e}*G*

_{up}.

#### 4) Water exchanges between aquifer and stream channel

The river package of MODFLOW simulates the water fluxes between stream channel and aquifers that are primarily controlled by the hydraulic conductivity and thickness of the bottom sediments and the head differences between streamflow and aquifer at the stream channel. The flow rate *Q _{g}* between the stream channel and aquifer is calculated by the difference in hydraulic heads in the stream and the adjacent aquifer using the following equation (McDonald and Harbaugh 1988):

where *C*_{riv} is the hydraulic conductance of the stream–aquifer interconnection, *h* is the head at the node in the cell underlying the stream reach, *H*_{riv} is the head in the stream channel. Obviously positive and negative values of *Q _{g}* represent groundwater recharge from river channel

*R*and baseflow discharge from aquifer

_{r}*R*, respectively.

_{g}### d. Routing of surface runoff and streamflow

Movement of surface runoff and streamflow is mainly regulated by the basin’s terrain and the channel’s characteristics. For runoff-generating storm events, surface runoff governs the peak discharge rate, time to peak, and volume of observed downstream outflows. Since the ground surface is relatively flat in a plains region, surface runoff generated in each grid cell flows into its nearby ditches and then concentrates in the stream channel. A simplified surface runoff routing method called a time lag approach, as described by the equation below, was used to represent the ground surface regulation:

where *Q _{s}*(

*t*) and

*Q*(

_{s}*t*− Δ

*t*) is the discharge at an outlet at time

*t*and

*t*− Δ

*t*, respectively,

*R*

_{s}is average value of surface runoff

*R*between time

_{s}*t*and

*t*− Δ

*t*, CS is coefficient of surface runoff concentration, and Lag is time lag.

Stream channels are divided into reaches for streamflow routing. Besides upstream inflow, each river reach may receive surface runoff and baseflow from the adjacent land segments. Stream discharge may also be influenced by water diversion, direct precipitation, and ET. Water balance analysis for each river reach provides an estimate of reach storage and thus its corresponding stage at each time step. The estimated storage and stage determine the discharge at each cross section, and streamflow is thus routed from upstream to downstream.

## 3. Model integration and implementation

The numerical modeling system is essentially an integration of the water flux algorithms described above, which are tightly coupled into the soil moisture model and the groundwater model. In practice, water fluxes are either inputs or outputs for the two models and thus they are numerically simulated through an iterative process. The physically based simulation of water fluxes among the atmosphere, stream channel, soil layers, and the aquifer is implemented using a finite-difference approach. A study site (watershed or field) is divided into rectangular or square cells, each of which is considered to be hydrologically and hydrogeologically uniform in land use, soil, and predominant depth of the water table.

Figure 2 schematically presents the water fluxes throughout the hydrological cycle and how they are interactively linked with reservoirs as represented in our modeling system. After precipitation is intercepted by vegetation a portion of it, *P _{d}* reaches the ground surface. Water may become surface runoff,

*R*, or infiltrate the soil and then further percolate into a groundwater aquifer,

_{s}*P*

_{rg}. Part of the water is returned to the atmosphere through evapotranspiration, ET. Driven by capillary force, groundwater may be pulled up to move into the bottom layer of the soil profile and then the moisture (

*G*

_{up}) may move farther upward as a result of a soil moisture gradient and finally return to the atmosphere through evapotranspiration. Human activities may introduce additional water fluxes into the hydrological cycle, namely, groundwater pumping,

*W*, stream water diversion for irrigation, and storing and percolation from artificial ponds.

_{g}Water fluxes interact with soil moisture dynamics, groundwater table fluctuations, and streamflow regimes. The modeling system simulates these interactions between each reservoir and their related water fluxes. Given the inputs of rainfall, potential ET, and land-cover and soil characteristics, the coupling of water flux algorithms with the soil moisture model can simulate the soil moisture dynamics and thus produce estimates of infiltration, surface runoff, groundwater recharge, and actual ET. For groundwater simulation, algorithms of the two vertical fluxes—groundwater recharge, *P*_{rg}, and loss, *G*_{up}—are inserted into the recharge (RCH) package and evapotranspiration (EVT) package of MODFLOW, respectively. The two water fluxes can be estimated by coupling the modified MODFLOW with the soil moisture model for iterative simulation based on the objective function of minimizing the errors between the observed and simulated values of streamflow and the groundwater table, as well as soil moisture contents if available. The horizontal flux of water exchanges between the aquifer and the stream channel is simulated directly by the river (RIV) package. Computational codes in FORTRAN and computer interfaces in Visual Basic have been developed to integrate all the algorithms into the modeling system.

The modeling system contains two types of parameters. The majority of parameters carry certain physical significance and they can be measured through experiments in the field or laboratory. These parameters primarily represent the soil and vegetation characteristics such as saturated hydraulic conductivity, water contents at wilting point, field capacity and saturation, and canopy capacity and reflectivity. Other empirical parameters such as the coefficient of surface runoff concentration, CS, time lag, storage capacity of small ponds, and hydraulic conductance, *C*_{riv}, need to be calibrated against the observed streamflow and the groundwater table. Finally, human influences on the water fluxes must be accounted for and irrigation is often the most important input into the soil water bucket or the output from an aquifer and stream channel. The amount of irrigation water can be estimated either by the difference between vegetation demand and the water available from rainfall and soil moisture storage or by using the records of water usage.

## 4. Model testing and application

### a. Study sites and model testing approach

The model was tested in two sites located in the Huaihe River basin where climate conditions vary from semihumid in the southeast to semiarid in the northwest (Fig. 1). The first site is the Linhuanji (LHJ) catchment with an area of 2 560 km^{2}. This catchment has meteorological and hydrological data covering a 10-yr period (1986–95), as well as spatial data of soil, vegetation, aquifer, and land characteristics. Therefore, it was chosen for model calibration and water budget analysis. Model validation was conducted in a runoff experimental field of 1.36 km^{2} surrounded by drainage ditches at the Wuduogou Hydrological Experimental Station (WHES). For this field, the available data include not only those required for model calibration against streamflow and the groundwater table for a 7-yr period (1994–2000), but also values of soil moisture content, as well as lysimeter-measured groundwater recharge and loss, which can be used to further evaluate the model’s accuracy.

In the LHJ catchment, average annual precipitation during 1986–95 is 713 mm and approximately 60%–70% of the precipitation occurs in the summer season from June to September. The annual potential evapotranspiration is 960 mm. Figure 3 shows the terrain location of groundwater observation wells and aquifer composition for this study basin. The ground surface elevations vary from 45 m above sea level in the north to 28 m above sea level in the south. The spatial distribution of the groundwater table demonstrates a similar pattern: 37–38 m above sea level in the north to 26–27 m above sea level in the south. The depth to groundwater table is as deep as 7–8 m in the north and decreases to only 2 m in the south. The groundwater table fluctuates about 1–2 m annually.

Data from the 10-yr period (1986–95) for the LHJ catchment includes daily precipitation recorded at 25 stations, pan evaporation from one station, groundwater table depths at 5-day intervals from 30 observation wells (Fig. 3b), and the daily streamflow discharge at the catchment outlet. The geological materials from aquifers in the northern part of LHJ, mostly sand and silty sand, came from the alluvial deposits of the abandoned Yellow River. In the southern part of the catchment, the aquifer sources are mostly silty sand, which is locally called “calcic concretionary black soil,” a main type of soil in this plain region of the Huaihe River basin (Fig. 3c). In the upper cultivated soils of unsaturated zone, sandy loam dominates in the north and sandy clay loam becomes more common in the south. Wheat, maize, and sorghum are the main crops in the region.

WHES has been a site for studies of the hydrological cycle in the Huaihe River plain region since the early 1950s. Regularly observed meteorological variables include precipitation, air temperature, soil temperature, humidity, wind speed, solar radiation, and pan evaporation. Other data collected at the experimental runoff field site include the water level in drainage ditches, the groundwater table, and soil moisture contents adjacent to the observation well at the depths of 5, 25, 55, and 90 cm below the ground’s surface. Observations indicate the average depth from the ground surface to the water table is 1.5 m. Additionally, 60 sets of lysimeters for soil columns of different depths were installed to measure groundwater recharge and loss. Soils at WHES represent a typical composition of the silty sand that is most common in the Huaihe River plain region. Wheat is the dominant crop at the experimental runoff field. Soil parameters for silty sand estimated in LHJ catchment were used for model validation in the field.

### b. Model calibration and validation

The LHJ catchment is discretized into 2356 grid units, each 1047 m long and 1048 m wide. The simulation time step for model calibration is one day. Results of the soil analysis by Song (2003) provide physical parameters of soil moisture dynamics, for example, volumetric water content at saturation *θ _{s}*, field capacity

*θ*, and wilting point

_{f}*θ*in the unsaturated zone. According to Cosby et al. (1984), hydraulic conductivity of sandy loam and sandy clay loam is 0.41 and 0.29 m day

_{w}^{−1}, respectively. In the saturated zone, hydraulic conductivity and the specific yield in the sandy region are 4.4 m day

^{−1}and 0.055 and are 2.8 m day

^{−1}and 0.045 in the silty sand region, respectively (Song 2003). Hydraulic conductivity of a riverbed deposit and channel cross-sectional characteristics determine the hydraulic conductance

*C*

_{riv}, ranging from 200 to 1500 m

^{2}day

^{−1}. Other parameters are calibrated against the observed stream discharge and the groundwater table using a trial and error method. The calibrated coefficient of the surface runoff concentration,

*CS*, and time lag is 0.35 and 1 day, respectively; the storage capacity of small ponds with an area of about 1% of the LHJ catchment is 2.5 mm.

The simulated and observed discharge for 1990–91 and 1994–95 examples are given in Fig. 4, showing that the two time series match each other reasonably well. The Nash–Sutcliff index for the entire calibration period (1986–95) is 0.79, and RMSE is 4.0 m^{3} s^{−1}. Model performance is also evaluated by comparing the simulated and observed groundwater table. As shown in Fig. 5, the groundwater table is simulated quite accurately for individual observation wells such as 13, 29, and 35. Averages of the simulated and observed water table for all of the 30 wells at each time step during the 10-yr calibration period are plotted in Fig. 6. Simulated and observed values of the water table at each well are averaged over the whole calibration period, and the correlation coefficient for the two sets of data is 0.81 (Fig. 7). These two figures further prove the model’s accuracy. Although the results of model calibration are generally satisfactory, it should be pointed out that some significant discrepancies between the observed and simulated values of discharge and the groundwater table do exist. We believe that these errors are unavoidable because the amount, timing, and location of human interferences on regional water fluxes and distribution, such as artificial pond storage and irrigation pumping, cannot be accurately estimated.

Model validation for the runoff experimental field in WHES also produces promising results. The correlation coefficient between the simulated and observed ditch stage is 0.80, and the correlation coefficient between the simulated and observed groundwater table is 0.85. The measurement of soil moisture contents, groundwater recharge, and loss at the field allows us to evaluate how accurately the model can simulate the water storage and fluxes within the soil–aquifer continuum. Figure 8 presents the observed and simulated soil moisture contents averaged for the four soil layers in each year from 1994 to 2000, indicating that the model simulation can capture the seasonal variations of soil moisture contents throughout the year. The annual mean recharge coefficient *P*_{rg}/*P* is the ratio of groundwater recharge to precipitation aggregated for one year. Corresponding to a certain groundwater depth, the coefficient of groundwater loss *G*_{up}/*E _{p}* is the ratio of the loss to potential ET averaged over the whole simulation period. As shown in Figs. 9 and 10, the simulated and observed values of these two coefficients match quite well in most cases except for the model warm-up period, which directly affects the result of the dry year, 1994. These results prove the accuracy and reliability of water flux simulation.

### c. Water budget analysis

After calibration and validation, the model can be used to simulate water fluxes and storage rates within the hydrological system to investigate the water transfer process and evaluate available water resources. For the LHJ catchment, simulation results of the annual water budget are given in Table 1. On average over the 10-yr period, the sum of precipitation and groundwater withdrawal that is primarily for agricultural utilization (column 2 + column 6 in Table 1) is 769.46 mm every year, and approximately 87% of this amount, that is, 669 mm (column 3), is lost through soil evaporation and vegetation transpiration. About 8.2% of precipitation, that is, 57.86 mm in column 9, becomes streamflow. In other words, the annual mean runoff coefficient equals 0.082 on average, and it varies from 0.02 in the dry year of 1994 to 0.24 in the wet year of 1991. About 15% of precipitation, that is, 103.60 mm (column 4), enters the groundwater aquifer, which means an annual mean recharge coefficient of 0.15, and it varies from 0.03 to 0.20. More than half of the recharge amount is lost either through evapotranspiration or as baseflow that discharges into the stream channel (column 5 + column 8), and the remaining amount stays in the aquifer. The total of the surface runoff and baseflow (column 7 + column 8) equals the simulated runoff (column 9). The yearly relative errors between the observed and simulated runoff fall within the range of 1.94%–22.84%. The majority of larger simulation errors occur primarily in the dry years because the artificial influences on the water balance become more significant. For the whole simulation period over the 10 yr, the relative error of simulated runoff is only 1.75%.

Water balance should be maintained for each storage reservoir of the hydrological system. Theoretically, input minus output for each reservoir should be equal to the storage change. Because of computational approximation, there might be a small difference between the two quantities, and this difference divided by the average of input and output is defined as balance error. Table 2 shows the 10-yr (1986–95) average water balance for unsaturated and saturated zones. For the unsaturated zone, the inputs include precipitation *P*, upward groundwater flux from the aquifer *G*_{up}, and irrigation water *I _{r}*, and the outputs include soil moisture loss through evapotranspiration ET, surface runoff

*R*, and groundwater recharge

_{s}*P*

_{rg}. As noted earlier, the soil profile is divided into only four layers for soil moisture simulation. This approximation has caused a small balance error of 0.37% as shown in Table 2. For the saturated zone, the inputs are the recharge from the bottom of the unsaturated zone

*P*

_{rg}, the streamflow recharge into aquifer

*R*, and the horizontal inflow into the aquifer from outside of the simulation domain

_{r}*B*

_{in}. Regarding the outputs from the aquifer, groundwater is lost through upward capillary flow into the soil layers (

*G*

_{up}), discharge into the stream as baseflow (

*R*), withdrawal (

_{g}*W*) for irrigation and other consumption (domestic and industrial water use), and the horizontal outflow (

_{g}*B*

_{out}) on the aquifer boundary. For the 10-yr simulation period, the difference between inputs and outputs into the saturated zone is exactly balanced by the storage change and therefore the error is zero in groundwater simulation using the MODFLOW.

## 5. Summary and conclusions

Estimation of water fluxes requires a modeling system that integrates algorithms for simulation of water movement among the atmosphere, stream channel, soil layers, and groundwater aquifer. For the Huaihe River basin where human activities, especially agriculture, have profound impacts on the hydrological cycle, such a tool is needed to understand the processes and evaluate the quantities of water being transferred from precipitation to evapotranspiration, surface runoff, soil moisture, groundwater, and streamflow. In this study, coupling a soil moisture model with the MODFLOW groundwater model has created a platform for us to integrate water flux algorithms into the numerical modeling system.

The modeling system is largely based on physical laws and employs a finite-difference method to numerically simulate water movement and fluxes in a horizontally discretized watershed or field. The majority of model parameters carry physical significance and can be determined by field and laboratory measurements or derived from watershed characteristics contained in GIS and remote sensing data. Several other empirical parameters need to be estimated by model calibration. The model has been successfully tested in two steps at two sites in the Huaihe River plain. First, it was calibrated in the Linhuanji catchment against observations of streamflow and the groundwater table in a 10-yr period (1986–95). Second, model validation against intensive measurements of soil moisture contents and groundwater fluxes was conducted at a very small runoff experimental field in the Wudouguo Hydrological Experimental Station. Calibration and validation results prove the accuracy and reliability of the modeling system.

This study has developed a useful tool for water resource managers to estimate water fluxes and evaluate water distribution in the Huaihe River plain under natural conditions and human perturbations. The physically based model offers promising applicability because it can be calibrated and customized with relative ease. The distributed-parameter design allows modelers to predict the changes of a regional water budget resulting from any scenarios of land use practices in different parts of the watershed. Water balance analysis based on model results will provide crucial information for the optimization of water allocation, irrigation scheduling, and other tasks in water resource assessment and planning.

## Acknowledgments

This research was supported by the Program for New Century Excellent Talents in University, China (NCET-04-0492), and partially supported by a grant from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project CUHK4627/05H).

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Xi Chen, State Key Laboratory of Hydrology–Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, China. Email: xichen@hhu.edu.cn

This article included in the The Global Energy and Water Cycle Experiment (GEWEX) special collection.