Integrated water system simulation by considering hydrological and biogeochemical processes : model development , with parameter sensitivity and autocalibration

Integrated water system modeling is a feasible approach to understanding severe water crises in the world and promoting the implementation of integrated river basin management. In this study, a classic hydrological model (the time variant gain model: TVGM) was extended to an integrated water system model by coupling multiple water-related processes in hydrology, biogeochemistry, water quality, and ecology, and considering the interference of human activities. A parameter analysis tool, which included sensitivity analysis, autocalibration and model performance evaluation, was developed to improve modeling efficiency. To demonstrate the model performances, the Shaying River catchment, which is the largest highly regulated and heavily polluted tributary of the Huai River basin in China, was selected as the case study area. The model performances were evaluated on the key water-related components including runoff, water quality, diffuse pollution load (or nonpoint sources) and crop yield. Results showed that our proposed model simulated most components reasonably well. The simulated daily runoff at most regulated and less-regulated stations matched well with the observations. The average correlation coefficient and Nash–Sutcliffe efficiency were 0.85 and 0.70, respectively. Both the simulated low and high flows at most stations were improved when the dam regulation was considered. The daily ammonium–nitrogen (NH4–N) concentration was also well captured with the average correlation coefficient of 0.67. Furthermore, the diffuse source load of NH4–N and the corn yield were reasonably simulated at the administrative region scale. This integrated water system model is expected to improve the simulation performances with extension to more model functionalities, and to provide a scientific basis for the implementation in integrated river basin managements.


Introduction
Severe water crises are global issues that have emerged as a consequence of the rapid development of social economy, and include flooding, water shortages, water pollution and ecological degradation.These crises have hindered the equitable development of regions by compromising the sustainability of vital water resources and ecosystems.It is impossible to address these crises within a single scientific discipline (e.g., hydrology, hydraulics, water quality or aquatic ecology) because of the complicated interactions among physical, chemical and ecological components of an aquatic ecosystem (Kindler, 2000;Paola et al., 2006).The paradigm of integrated river basin management may be a sensible solution at basin scale by focusing on the coordinated management of water resources in terms of social economy, water quality and ecosystems.Integrated water system models have been popular since the last decade due to the rapid de- The hydrological cycle has been known as a critical linkage among other water-related processes (e.g., physical, biogeochemical and ecological processes) and energy fluxes at the basin scale (Burt and Pinay, 2005).For example, physiological and ecological processes of vegetation affect evapotranspiration, soil moisture distribution and nutrient movement.In the meantime, soil moisture and nutrient constrain the vegetation growth.Overland flow is a carrier of pollutants to water bodies.Therefore, all the processes should be considered simultaneously to capture the interactions and feedbacks between individual cycles.Multidisciplinary research provides an effective way to enable breakthroughs in the integrated water system modeling by integrating the theories in water-related sciences (e.g., accumulated temperature law for phenological development, Darcy's law for groundwater flow, Saint-Venant equation for flow routing, balance equation for mass and momentum, Richards' equation for unsaturated zone, Horton theory for infiltration, Penman-Monteith equation for evapotranspiration).Abundant open data sources further support the implementation of an integrated water system model, e.g., high-resolution spatial information data, chemical and isotopic data from field experiments (Singh and Woolhiser, 2002;Kirchner, 2006).
Several models have been developed since the 1980s (Di Toro et al., 1983;Brown and Barnwell, 1987;Johnsson et al., 1987;Hamrick, 1992;Li et al., 1992;Abrahamsen and Hansen, 2000;Tattari et al., 2001;Singh and Woolhiser, 2002).Owing to the complexity of the integrated water system and the scale conflicts between different processes, most existing models focus on only one or two major water-related processes, and can be categorized into three major classes.(1) Hydrological models emphasize the rainfall-runoff relationship and link with some dominating water quality and biogeochemical processes.These models generally show satisfactory performances in simulating the hydrological processes.Some widely accepted models are TOPMODEL (Beven and Kirkby, 1979), SHE (Abbott et al., 1986), HSPF (Bicknell et al., 1993), VIC (Liang et al., 1994), ANSWERS (Bouraoui and Dillaha, 1996), HBV-N (Arheimer andBrandt, 1998, 2000), HYPE (Lindström et al., 2010) and its improved version S-HYPE (Strömqvist et al., 2012).( 2) Water quality models focus on the migration and transformation processes of pollutants in water bodies.These models can simulate the water quality variables at high spatial and temporal resolutions in river networks by adopting multi-dimensional dynamic equations.However, they have difficulties in simulating the overland processes of water and pollutants.Typical models include WASP (Di Toro et al., 1983), QUAL2E (Brown and Barnwell, 1987) and EFDC (Hamrick, 1992).(3) Biogeochemical models have advantages in simulating the physiological and ecological processes of vegetation and the vertical movements of nutrients and water in soil layers at the field or experimental catchment scales.However, these models lack accurate hydrological features (Deng et al., 2011) and are hard to simulate the movements of water, nutrients and their losses along flow pathways in the basin.Some biogeochemical models are SOILN (Johnsson et al., 1987), EPIC (Sharpley and Williams, 1990), DNDC (Li et al., 1992), Daisy (Abrahamsen and Hansen, 2000) and ICECREAM (Tattari et al., 2001).Overall, most models usually achieve good performances on their oriented processes and only approximate the results for other processes outside of the model's focus in the integrated river basin management.An important scientific question is "does including these extra processes in an integrated manner improve model results compared to models that are focused only on one component?" SWAT is an integrated water system model that can simulate most water-related processes over a long period at large scales (Arnold et al., 1998).However, not all water-related processes can be well captured in practice because of the inaccurate descriptions of some processes, such as daily extreme flow events (Borah and Bera, 2004), soil nitrogen and carbon (Gassman et al., 2007) and regulation rules of dams or sluices in regulated basins (Zhang et al., 2013).Particularly, the simulation methods of surface runoff yield in SWAT have been questioned, e.g., the general applicability of the curve number (Rallison and Miller, 1981) and the scale limitations of the Green-Ampt infiltration model (King et al., 1999).Furthermore, SWAT has difficulties in accurately capturing the complicated dynamic processes of soil nitrogen and carbon by comparing with other biochemical models (Gassman et al., 2007).Several modified versions have been developed, such as SWIM (Krysanova et al., 1998) and SWAT-N (Pohlert et al., 2006).
In this study, we tended to develop an integrated water system model based on a hydrological model.The time variant gain model (TVGM) proposed by Xia (1991) is a lumped hydrological model based on the rainfall and runoff observations from many basins with different scales all over the world.In the TVGM, the rainfall-runoff relationship is considered to be nonlinear because the surface runoff coefficient varies over time and is significantly affected by antecedent soil moisture.The TVGM has a strong mathematical basis because this nonlinear relationship is transformed into a complex Volterra nonlinear formulation.Wang et al. (2002) extended the TVGM to the distributed time variant gain model (DTVGM) by taking advantage of better computing facilities and available data sources.Currently, the DTVGM performs well in many basins with different scales and climate zones to investigate the effect of human activities and climate change on runoff (Xia et al., 2005;Wang et al., 2009).
In the model development, we would like to produce reasonable simulations simultaneously in both hydrological and water quality processes and to include more water-related processes such as soil biogeochemistry and crop growth for better understandings of the complicated water-related processes and their interactions in the real basins.Our proposed model was built by extending the DTVGM through coupling of the detailed interactions and linkages among hydrological, water quality, soil biogeochemical and ecological processes, as well as considering the prevalent regulations of water projects (dams and sluices) at the basin scale.In order for readers to use the proposed model easily, a parameter analysis module, which included popular objective functions, autocalibration approaches and summary statistics, was also developed.To demonstrate the model performances, we simulated several key water-related components including flow regimes, diffuse source (or nonpoint source) pools of nutrients, water quality variables in water bodies and crop yield in a highly regulated and heavily polluted catchment (Shaying River catchment) in China.

Model framework
Our proposed model includes eight major modules, namely the hydrological cycle module (HCM), soil biochemical module (SBM), crop growth module (CGM), soil erosion module (SEM), overland water quality module (OQM), water quality module of water bodies (WQM) and dam regulation module (DRM).The parameter analysis tool (PAT) is also designed for model calibration.The model structure is shown in Fig. 1.More detailed descriptions of each module and its interactions with other modules are given in Sects.2.1.1 to 2.1.5.The main equations of each module are deferred to the Appendix and Supplement for readers who are interested in the mathematical details.
Our model is based on the hypothesis that the cycles of water and nutrients (N, P and C) are inseparable and act as the critical linkages among all the modules.It takes full advan-tage of the existing models, i.e., the powerful interconnections of the hydrological models with other processes at the spatial scale, the elaborative descriptions of the ecological models on nutrient vertical movement in soil layers, and the elaborative descriptions of the water quality models on nutrient movements along river networks.First, several key components, simulated by the hydrological cycle module (HCM) (e.g., evapotranspiration, soil moisture and flow), are treated as critical linkages in all the modules (Sect.2.1.1).Second, the soil biochemical processes determine the nutrient loads absorbed in the crop growth process (CGM) and migrated into water bodies as the diffuse pollution source (OQM and WQM).The accurate descriptions of soil biochemical processes are helpful in improving the simulation of diffuse source processes in responding to agricultural management (Sect.2.1.2).Third, the hydrological cycle module (HCM) provides a function for describing the connections between spatial calculation units to simulate the overland and instream movements of water and nutrients at the basin scale (Sects.2.1.1 and 2.1.3).

Hydrological cycle module (HCM)
Surface runoff calculation is the core of hydrological simulation.The TVGM is adopted to calculate the surface runoff yields for different land-use/cover areas, such as forest, grassland, water body, urban area, unused land, paddy land and dryland agriculture.The potential evapotranspiration is calculated using the Hargreaves method (Hargreaves and Samani, 1982) because only the available daily maximum and minimum temperatures are used.The actual plant transpiration is expressed as a function of potential evapotranspiration and leaf area index, whereas soil evaporation is expressed as a function of potential evapotranspiration and  surface soil residues (Neitsch et al., 2011).The yields of interflow and baseflow have linear relationships with the soil moisture in the upper and lower layers, respectively (Wang et al., 2009).The infiltration from the upper to lower soil layers is calculated using the storage routing method (Neitsch et al., 2011).The Muskingum method or kinetic wave equation is used for river flow routing.
Figure 2 shows that the shallow soil moisture from the hydrological cycle module is a major factor that connects the crop growth module (to control crop growth) and the soil biochemical module (to control the vertical migration and reaction of nutrients in the soil layers).Plant transpiration is also linked to the soil biochemical module (to drive the vertical migration of nutrients in the soil layers).The surface runoff is linked to the soil erosion module, while the overland flow (surface runoff, interflow and baseflow) is connected to the overland water quality module (to drive the movements of nutrients and sediment along flow pathways) and the water quality module of water bodies (rivers and lakes) for runoff routing.Moreover, the hydrological cycle module provides the inflows for individual dams or sluices in the dam regulation module.

Modules for ecological processes
The ecological processes are described in the soil biochemical module and the crop growth module.The crop growth and soil biochemical processes directly affect the soil moisture, evapotranspiration, nutrient transformation and loss from the soil layers.Therefore, our model incorporates the water cycle, nutrient cycle, crop growth and their key linkages.

Soil biochemical module (SBM)
The soil biochemical module simulates the key processes of carbon (C), nitrogen (N) and phosphorus (P) dynamics in the soil layers, including decomposition, mineralization, immobilization, nitrification, denitrification, leaching and plant uptake.Different forms of N and P outputted from the soil bio- chemical module are connected to the crop growth module as the nutrient constraints of crop growth and to the overland water quality module as the main diffuse sources to water bodies (Fig. 3a).
Soil C and N cycle.The sub-models of daily step decomposition and denitrification in DNDC (Li et al., 1992) are adopted to simulate the soil biogeochemical processes of C and N at the field scale.The decomposition and other oxidation processes are the dominant microbial processes in the aerobic condition.The three conceptual organic C pools are the decomposable residue C pool, microbial biomass C pool and stable C pool.The decomposition of each C pool is treated as the first-order decay process with the individual decomposition rates constrained by the soil temperature and moisture, clay content and C : N ratio.The major simulated processes of decomposition under aerobic conditions are mineralization, immobilization, ammonia (NH 3 ) volatilization and nitrification.The mineralization and immobilization of mineral N (NH + 4 and NO − 3 ) are determined by the flow rates of soil organic carbon (SOC) pools.NH 3 volatilization is controlled by the NH + 4 concentration, clay content, pH, soil moisture and temperature.NH + 4 is oxidized to NO − 3 during nitrification and nitrous oxide (N 2 O) is emitted into the air during the nitrification.Denitrification occurs under the anaerobic condition, which is controlled by soil moisture, temperature, pH and dissolved SOC content.The detailed descriptions are given in Appendix B and Li et al. (1992).
Soil P cycle.The major processes of the soil P cycle are simulated according to the study of Horst et al. (2001).Six P pools are considered including three organic pools (stable and active pools for plant uptake, a fresh pool associated with plant residue) and three mineral pools (dissolved mineral, stable and active pools).The involved processes are the P release, mineralization and decomposition from fertilizer, manure, residue, microbial biomass, humic substances and the sorption by plant uptake (Horst et al., 2001;Neitsch et al., 2011).
The soil profile is divided into three layers, namely, surface (0-10 cm) and user-defined upper and lower layers, all of which are consistent with the soil layers of the hydrological cycle module to smoothly exchange the values through the linkages (e.g., soil moisture) among different modules.

Crop growth module (CGM)
The crop growth module is developed based on the EPIC crop growth model (Hamrick, 1992).It simulates total dry matter, leaf area index, root depth and density distribution, harvest index, nutrient uptake, and so on (Williams et al., 1989;Sharpley and Williams, 1990).The crop respiration and photosynthesis drive the vertical movements of water and nutrients.The output of the leaf area index is a main factor connecting the hydrological cycle module (to control the transpiration), and the crop residue left in the fields is a main source of organic nutrients (C, N and P) connecting to the soil biochemical module for soil biochemical processes, to the overland water quality module and to the soil erosion module as one of the five constraint factors (Fig. 3b).

Modules for water quality processes
The water quality processes focus on the migration and transformation of water quality variables (e.g., sediment, different forms of nutrients, biochemical oxygen demand, BOD, and chemical oxygen demand, COD) along the flow pathways in the land surface and river network.The main modules are the soil erosion module for the sediment yield, the overland water quality module for the migration of overland diffuse source to water bodies and the water quality module for the migration and transformation of point and diffuse pollution sources in water bodies.

Soil erosion module (SEM)
The soil erosion by precipitation is estimated using the improved USLE equation (Onstad and Foster, 1975) based on runoff yields outputted from the hydrological cycle module and crop management factor outputted from the crop growth module.The soil erosion module simulates the sediment load for the overland water quality module to provide the carrier for the migration of insoluble organic matter along overland transport paths and water bodies (Fig. 4a).

Overland water quality module (OQM)
This module simulates the overland losses and migration loads of diffuse source pollutants (e.g., sediment, insoluble and dissolved nutrients, BOD and COD) (Fig. 4b).The main diffuse sources include the nutrient loss from the soil layers and urban areas, and the farm manure from livestock in rural areas.The nutrient loss from the soil layers, as the primary diffuse source in most catchments, is determined by the overland flow and sediment yield (Williams et al., 1989), and the other sources are estimated using the export coefficient method (Johnes, 1996).The overland migration processes contain the dissolved pollutant migration with overland flow and the insoluble pollutant migration with sediment.All the processes occur along the overland transport paths.

Water quality module of water bodies (WQM)
This module simulates the transformation and migration of water quality variables in different types of water bodies (in-stream and water impounding) (Fig. 4c).The simulated variables include water temperature, dissolved oxygen (DO), sediment, different forms of nutrients (N and P), Two modules are designed for the different types of water bodies, i.e., the in-stream water quality module and the water quality module for water impounding (reservoir or lake).The enhanced stream water quality model (QUAL-2E) (Brown and Barnwell, 1987) is adopted to simulate the longitudinal movement and transformation of water quality variables in the in-streams.The model is solved at the sub-basin scale rather than at the fine grid scale in order to maintain spatial consistency with the hydrological cycle module.The water quality outputs provide the water quality boundary of dams or sluices in the dam regulation module.The water quality module for water impounding assumes that water body is at the steady state and focuses on the vertical interaction of water quality processes.The main processes include water quality degradation and settlement, sediment resuspension and decay.

Dam regulation module (DRM)
Dams and sluices highly alter flow regimes and associated water quality processes in most river networks.Thus, the dam and sluice regulation should be considered in the water system models.The dam regulation module provides the regulated boundaries (e.g., water storage and outflow) to the hydrological cycle module for flow routing and to the water quality module of water bodies for pollutant migration.
Given that different types of dams and sluices are likely to show completely different regulation behaviors, we try to reproduce their common functionalities for either the flood control or water supply in this module.Three methods are proposed to calculate the water storage and outflow of dams or sluices, namely, the measured outflow, controlled outflow with target water storage and the relationship between outflow and water storage volume.The first method requires users to provide the measured outflow series during the simulation period.The second method simplifies the regulation rules of dams or sluices for long-term analysis based on the assumption that water is stored according to the usable water level during the non-flooding season and the flood control level during the flooding season, and the surplus water is discharged.This method requires the characteristic parameters of dam or sluice including water storage capacities of dead, usable, flood control and maximum flood levels and the corresponding water surface areas.The third method is based on the relationships among water level, water surface area, storage volume and outflow according to the designed dam data or long-term observed data (Zhang et al., 2013) (Appendix C).

Parameter analysis tool (PAT)
In our model, 66 lumped and 94 distributed parameters involve the hydrological, ecological and water quality processes.The distributed parameters are divided into 37 overland parameters, 17 stream parameters and 40 parameters of water projects (only for the sub-basin with reservoir or sluice) according to their spatial distribution.These parameter values are determined by the properties of overland landscape and soil, stream patterns and water projects, respectively.Different spatial calculation units share many common parameter values if their properties are the same.
Owing to a large number of parameters, it is hard to find optimal parameter values by manual tuning.The limited number of observed processes causes equifinality in the model calibration (Beven, 2006).Therefore, the parameter sensitivity analysis and calibration are important steps to alleviate equifinality in the applications of highly parameterized models, particularly for integrated water system models (Mantovan and Todini, 2006;Mantovan et al., 2007;McDonnell et al., 2007).The PAT is designed for parameter sensitivity analysis, autocalibration and model performance evaluation (Fig. 5).
To evaluate model performance, five traditionally used criteria are included in the PAT, i.e., bias (bias), relative error (re), root mean square error (RMSE), correlation coefficient (r) and Nash-Sutcliffe efficiency (NS defined by Nash and Sutcliffe, 1970).The detailed definitions of these criteria are given in Appendix D. Furthermore, flow duration curve and cumulative distribution function are also provided for capturing multiple signatures of calibrated processes.More criteria can also be proposed by the users.The objective function(s) to calibrate the model can be formed by single or multiple criteria or their function (such as weighted average).
The parameter analysis algorithms in the PAT include the parameter sensitivity method (Latin hypercube one factor at a time: LH-OAT) ( van Griensven et al., 2006), the single objective auto-optimization methods such as particle swarm optimization (PSO) (Kennedy, 2010), the genetic algorithm (GA) (Goldberg, 1989) and shuffled complex evolution (SCE-UA) (Duan et al., 1994), as well as the multi-objective autooptimization methods such as the weighted sum method and nondominated sorting genetic algorithm II (NSGA-II) (Deb et al., 2002).The method can be selected on the basis of the specific requirements of users.
In order to obtain the optimal parameter values, the following treatments are adopted in the PAT.First, the prior ranges of all the parameter values or their prior distributions (i.e., uniform or normal) are preset by referring to the literatures or similar basins.The constraints on parameters are also considered in both parameter sensitive analysis and autocalibration.In the hydrological cycle module, the constraints on soil moisture parameters are W m (minimum moisture) < W w (moisture at permanent wilting point) < W fc (field capacity) < W sat (saturated moisture capacity).The basic surface runoff coefficient (g 1 ) for different land uses/covers is set in ascending order (water body, paddy land, urban area, forest, dryland agriculture, unused land and grassland).The interflow yield coefficient (K ss ) is greater than the baseflow coefficient (K bs ).In the water quality module of water bodies, the settling rates of water quality variables (K set ) in the water impounding are greater than the resuspension rates (K scu ) and the settling rates (R set ) in channels.Second, the sensitive parameters are determined to reduce the parameter dimensions by sensitivity analysis.Third, the selected sensitive parameters are calibrated by the auto-optimization method, while the insensitive parameters remain as their default values that are given by referring to the literatures or other models (e.g., SWAT, EPIC and DNDC) in the same/similar basins.
The PAT connects with other modules through the parameter values that are used to simulate the processes of other modules and evaluate the objective functions in sensitivity analysis and autocalibration.Depending on the algorithm used, the parameter values are (randomly) sampled from the multi-dimensional parameter spaces to drive our model, and the objective function value of each parameter set is then obtained.For the parameter sensitivity analysis, the sensitivity index of each parameter set is evaluated by comparing the variation of the objective function value along with the change in parameter value.For the parameter autocalibration, the good parameter sets are kept or updated by the autooptimization method until the convergence or the maximum number of iterations is achieved.

Multi-scale solution
The spatial heterogeneities of basin attributes and the different timescales used in individual processes cause inconsistent spatial and temporal scales in model integration (Sivapalan and Kalma, 1995;Singh and Woolhiser, 2002).For the spatial scale, three levels of spatial calculation units are designed, namely, sub-basin, land-use/cover and crop from largest to smallest.These units are defined as the minimum polygons with similar hydrological properties, land uses/covers and agriculture crop cultivation patterns, respectively.The sub-basins are defined on the basis of a digital elevation model (DEM), the positions of gauges and water projects, and are used in the hydrological cycle module (e.g., flow routing in both land and in-stream), overland water quality module, water quality module of water bodies and dam regulation module.Seven specific land-use/cover units of each sub-basin are partitioned by the land-use/cover classification (i.e., forest, grassland, water, urban, unused land, paddy land and dryland agriculture) and are used in the hydrological cycle module (e.g., water yield, infiltration, interception and evapotranspiration) and the soil erosion module.Moreover, several specific land-use/cover units (paddy land, dryland agriculture, forest and grassland), where agricultural activities usually occur, are divided further into the crop units for the detailed analysis of the impact of agricultural management on water and nutrient cycles.In the current version of our model, these four land-use/cover units are divided into 10 specific categories of crop units as fallow for all these land-use/cover units, grass for grassland unit, fruit tree and non-economic tree for forest unit, early rice and late rice for paddy unit, spring wheat, winter wheat, corn and mixed dry crop for dryland agriculture unit.The crop unit of a specific land-use/cover pattern varies depending on crop cultivation structure and timing.The related modules are the soil biochemical module and the crop growth module.All of the outputs of the crop unit are summarized at the land-use/cover scale or sub-basin scale based on the area percentages in different crop units.
For the temporal scale, it is practical to use a daily time step, as this is consistent with the underlying rainfall-runoff module and the data availability.The sub-daily scale may improve the performance in some modules (e.g., SEM and WQM).However, most observations (e.g., climate data sets, soil nutrient availability and water quality concentrations) are at the daily scale, leading to potential uncertainties or instabilities to disaggregate the observations into a sub-daily scale.Linear or nonlinear aggregation functions are used to transform different timescales to daily scale (Vinogradov et al., 2011), such as exponential functions for flow infiltration and overland flow routing processes in the hydrological cycle module, for soil erosion processes in the soil erosion module (Eqs.A5, A6 and S32 in Appendix A and the Supplement), and accumulation functions for the crop growth process in the crop growth module (Eq.S7 in the Supplement).

Basic data sets and spatial delineation
The indispensable data sets for model setup are GIS data, daily meteorological data series, social and economic data series and dam attribute data.Several monitoring data series are needed for model calibration, such as runoff and water quality series in river sections, soil moisture and crop yield at the field scale.Table 1 shows all of the detailed data sets and their usages.
The hydrological toolset of the Arc GIS platform is used to delineate all the spatial calculation units based on a DEM and land-use/cover data.The sub-basin attributes (e.g., location, evaluation, area, land surface slope and slope length, landuse/cover areas) and flow routing relationship between subbasins are obtained during this procedure.

Study area and model testing
In this study, our model was applied to a highly regulated and heavily polluted catchment (the Shaying River catchment) in China.The simulated water-related components contained daily runoff and water quality concentrations at river sections, spatial patterns of diffuse source pollution load and crop yield at sub-basin scale.

Study area
The Shaying River catchment (112 • 45 -113 • 15 E, 34 • 20 -34 • 34 N), which is the largest sub-basin of the Huai River basin in China, is selected as the study area (Fig. 6a).The drainage area is 36 651 km 2 , with a mainstream of 620 km.The average annual population (2003-2008) (Fig. 6b) is 32.42 million, with a rural population of 23.70 million.The average annual stocks include 8.30 million big animals (cattle, pigs and sheep) and 178.42 million poultry (Fig. 6c).The average annual use of chemical fertilizer is 1.55 million ton (N: 38-51 %; P: 16-25 %; and others: 23-47 %) (Fig. 6d).The catchment is located in the typical warm temperate and Hydrol.Earth Syst.Sci., 20, 529-553, 2016 www.hydrol-earth-syst-sci.net/20/529/2016/ semi-humid continental climate zone.The annual average temperature and rainfall are 14-16 • C and 769.5 mm, respectively.The Shaying River is the most seriously polluted tributary, with a pollutant load contribution of over 40 % in the whole Huai River, and is usually known as the water environment barometer of the Huai River mainstream.To reduce flood or drought disasters, 24 reservoirs and 13 sluices, whose regulation capacities are over 50 % of the total annual runoff, have been constructed, and fragmented the river into several impounding pools.

Model setup
All data sets for model setup and calibration were collected from the government bureaus, official books and scientific references.The detailed descriptions were presented in Tables S2 and S3 of the Supplement.The resolutions of GIS and weather input data were quite satisfactory for the model application.However, most data on water quality, ecology and agricultural management were at monthly or annual temporal scale.The data for economy, agricultural management and diffuse source load were collected from individual administrative regions.Both the temporal and spatial scales were larger than the required daily scale or spatial calculation units (sub-basin, land-use/cover and crop).In these cases, the data values were uniformly distributed to the required temporal and/or spatial scales, such as the input of point sources, and social and economic data.The Shaying River catchment was divided into 46 subbasins.According to the land-use/cover classification standard of China (CNS, 2007), the main land-use/cover types were dryland agriculture (84.04 %), forest (7.66 %), urban (3.27 %), grassland (2.68 %), water (1.43 %), paddy land (0.91 %) and unused land (0.01 %).The soil input parameters (the contents of sand, clay and organic matter) were calculated based on the percentage of soil types in each subbasin.The main crops were early rice and late rice in the paddy land, and winter wheat and corn in the dryland agriculture.The main agricultural management schemes (fertilize, plant, harvest and kill) were summarized by field investigation in the studies of Wang et al. (2008) and Zhai et al. (2014) (Table S3).Crop rotations and management schemes were considered in the model by setting the start time, the duration of management and the fertilizer amounts.Two fertilizations (base and additional fertilization) were considered in the model during the complete growth cycle of a certain crop.The areas of sub-basin, land-use/cover and crop units ranged from 46.48 to 3771.15 km 2 , from 0.04 to 2762.5 km 2 , and from 3.73 to 2762.5 km 2 , respectively.
The daily precipitation series from 2003 to 2008 at 65 stations were interpolated to each sub-basin using the inverse distance weighting method, while the daily temperature series at six stations were interpolated using the nearestneighbor interpolation method.The social and economic data (e.g., population and livestock in the rural area, chemical fer-tilizer amounts) were calculated for each sub-basin based on the area percentage.
Moreover, 5 reservoirs, 12 sluices and over 200 wastewater discharge outlets were considered according to their geographical positions.The farm manure from rural living and livestock farming was considered as a diffuse source owing to its scattered characteristics and the deficient sewage treatment facilities in the rural areas.

Model evaluation
The observation series of daily runoff and NH 4 -N concentration were used to calibrate the model parameters.There were five regulated stations (Luohe, Zhoukou, Huaidian, Fuyang and Yingshang) and one less-regulated station (Shenqiu), which is the downstream station situated far from water projects.Moreover, given that the observed yields of diffuse pollutant loads and crops were hard to collect for the whole catchment, only the statistical results from official reports or statistical yearbooks (Wang, 2011;Henan Statistical Yearbooks, 2003, 2004and 2005) were collected to validate the model performances.
We selected LH-OAT for parameter sensitivity analysis and SCE-UA for parameter calibration in the PAT.To reduce the dimensions of the calibration problem, we restricted SCE-UA to calibrate only the sensitive parameters defined by LH-OAT, whereas the rest of the parameters remained constants.The selected evaluation indices of model performance were bias, r and NS.However, NS was sensitive to the extreme value, outlier and number of the data points, and was not commonly used in environmental sciences (Ritter and Muñoz-Carpena, 2013).Thus NS was not used to evaluate the NH 4 -N concentration simulation.
The model calibration was conducted by the following steps.Hydrological parameters were calibrated first against the observed runoff series at each station from upstream to downstream, and then water quality parameters against the observed NH 4 -N concentration series.The calibration and validation periods were from 2003 to 2005 and from 2006 to 2008, respectively.The weighted sum method was usually used to comprehensively handle multi-objectives (Efstratiadis and Koutsoyiannis, 2010).In this study, singleobjective functions were formed by equally weighting the evaluation indices as (f runoff and f NH4−N ) because the case study was only a demonstration of the model performance.
Moreover, the effect of dam regulation was considered because of the high regulation in most rivers.The dam and sluice regulation usually altered the intra-annual distribution of flow events, such as flattening high flow and increasing low flow.The simulation performances of high and low flows were separately evaluated and the effectiveness of the DRM was tested by comparing the simulation with and with- out the consideration of dam regulation.The high and low flows were determined by the cumulative distribution function (CDF).A threshold of 50 % was used for easy presentation; i.e., the flow was treated as high flow (or low flow) if its percentile was greater than (or smaller than) the threshold.

Parameter sensitivity analysis
Nine sensitive parameters were detected for runoff simulation by LH-OAT (Table 2), including soil-related parameters W fc (field capacity), W sat (saturated moisture capacity), K r (interflow yield coefficient) and K sat (steady-state infiltration rate); TVGM parameters g 1 (basic surface runoff coefficient) and g 2 (influence coefficient of soil moisture); baseflow parameters K g (baseflow yield coefficient) and T g (delay time for aquifer recharge); and evapotranspiration parameter K ET (adjusted factor of actual evapotranspiration).All of these parameters controlled the main hydrological processes in which soil water and evapotranspiration processes were distinctly important and explained 54.3 and 23.2 % of the runoff variation, respectively.For NH 4 -N concentration simulation, over 90 % of observed NH 4 -N concentration variations were explained by 14 sensitive parameters that were categorized into hydrological (59.28 % of variation), NH 4 -N (20.65 % of variation) and COD (12.34 % of variation) related parameters.The main explanation was that hydrological processes provided the hydrological boundaries that affected the diffuse source load into rivers and the degradation and settlement processes of NH 4 -N in water bodies NH 4 -N concentration was further influenced by the settlement and biological oxidation.More-over, it was a competitive relationship between COD and NH 4 -N to consume DO of water bodies in a certain limited level (Brown and Barnwell, 1987).

Hydrological simulation
The runoff simulations fitted the observations well at all the stations (Fig. 7 and Table 3).The biases were very close to 0.0 at all the regulated stations except Zhoukou with an underestimation (bias: 0.24 for calibration and 0.41 for validation) and Luohe with an overestimation (bias: −0.52 for validation).The obvious biases were caused by the average objective function of all three evaluations rather than the bias only.The r values ranged from 0.75 (Luohe for validation) to 0.92 (Yingshang for calibration) with the average value of 0.85, whereas the NS values ranged from 0.51 (Luohe for validation) to 0.84 (Yingshang for calibration) with the average value of 0.70.The results of the regulated stations were a little worse than those of the less-regulated station (Shenqiu) owing to the regulation.
By comparing the simulations with the observations from 2003 to 2008, we saw that the high and low flows were always overestimated if the model did not consider the regulations (Fig. 8).Except for the high flows at Zhoukou, both high and low flows at all the stations were simulated well when the dam and sluice regulation was considered (Table 4).The best fitting was at Fuyang, particularly for the high flow simulation (bias = 0.10, r = 0.89 and NS = 0.78).From unregulation to regulation settings, the improvements measured by f runoff ranged from −0.08 (Zhoukou) to −0.29 (Huaidian) for high flow simulations, from −0.05 (Zhoukou) to −0.31 (Huaidian) for average flow simulations, and from  flow simulations were very obvious.However, their performances still needed to be improved further, particularly for the underestimation at Zhoukou and Huaidian.The possible reasons were as follows.On the one hand, the applied evaluation indices (r and NS) were known to emphasize the high flow simulation rather than the low flow simulation (Pushpalatha et al., 2012), and the objective of autocalibration was to obtain the optimal solution for the average of three evalu-ation indices rather than the bias only.The slight sacrifice of bias improved the overall simulation performance evaluated by all three indices.One the other hand, the dam regulation module still could not fully capture the low flows.Furthermore, the model performances on monthly flows were even better, particularly for r and NS.The r values ranged from 0.87 (Luohe for both calibration and validation) to 0.95 (Fuyang for calibration) with the average value of Table 4.The runoff simulation results at regulated stations with and without the dam regulation considered.Range means the difference of the objective function value between regulations considered and not considered.If the range value is less than 0.0, then the simulation with regulation is better than that without regulation.Otherwise, the simulation without regulation is better.

Stations
Regulated capacity (%) Flow event 0.92, whereas the NS values ranged from 0.67 (Luohe for validation) to 0.94 (Shenqiu for validation) with the average value of 0.80.Compared with the existing results at the same stations by SWAT (Zhang et al., 2013), the flow simulations at the downstream stations were improved, although they became a little worse at the upstream stations (Luohe and Zhoukou for calibration).In particular, the total water volume and agreements with the observations (i.e., bias and NS) were well captured.

Water quality simulation
The simulated concentrations of NH 4 -N matched well with the observations according to the evaluation standard recommend by Moriasi et al. (2007) (Fig. 9 and Table 5).The r values were over 0.60 for all the stations except Zhoukou (0.56 for validation), Yingshang (0.49 for validation) and Shenqiu (0.41 for validation), and the average value was 0.67.The biases were considered to be "acceptable" with a range from −0.27 (Fuyang for validation) to 0.29 (Zhoukou for calibration).The best simulation was at Luohe station.The obvious discrepancies between the simulations and observations often appeared in the period from January to May because of the poor simulation performances on the low flows.Although the biases changed markedly from calibration to validation at Fuyang and Yingshang stations, the model performances were still acceptable.The possible explanation was that the biases for corresponding runoff simulations at these two stations also changed.
Compared with the results without the consideration of regulation, the simulation results were obviously improved when the regulation was considered, except those at Fuyang station in the calibration period.The decreases in the f NH4−N value ranged from 0.10 (Huaidian for calibration) to 0.49 (Zhoukou for validation), although there was a slight increase at Fuyang for the calibration (0.02).Therefore, it was concluded that the consideration of dam and sluice regulation played an important role in the water quality simulation.In the upper stream of the Shaying River, the flow was small and the NH 4 -N concentration decreased obviously because of the degradation and settlement of large water storage.In the downstream of the Shaying River, the NH 4 -N concentration increased because of the pollutant accumulation and the decreasing flow from dams and sluices owing to the regulation (Zhang et al., 2010).Therefore, the simulated concentrations without regulation were usually overestimated or higher than the simulation with regulation at the upstream stations (Luohe and Zhoukou).However, the concentrations were underestimated at the downstream stations (Huaidian, Fuyang and Yingshang).The largest differences between the simulations with and without the consideration of regulation appeared at Zhoukou.
The spatial pattern of average annual load of diffuse source NH 4 -N was shown in Fig. 10a.The estimated annual yield rates ranged from 0.048 to 11.00 t km −2 year −1 with the average value of 0.73 t km −2 year −1 .The yield in each administrative region was summarized from the results of each subbasin according to the area percentage of sub-basins in each administrative region.Compared with the statistical load of each administrative region based on the soil erosion, land use/cover and fertilizer amount in the official report (Wang, 2011) region was 21.31 % when the two regions with the biggest biases (Fuyang and Pingdingshan) were excluded as outliers.
The high load regions were in the middle of the Pingdingshan, Xuchang, Zhengzhou, Fuyang and Zhoukou regions.
The spatial pattern was significantly correlated with the distribution of paddy area (r = 0.506, p < 0.001) and rice yield (r = 0.799, p < 0.001) (Fig. 10b and c).The fertilizer losses in the paddy areas might be the primary contributor to the diffuse source NH 4 -N load because the average nitrogen loss coefficient in China was just 30-70 % in the paddy areas, which was higher than that in the dryland agriculture (20-50 %) (Zhu, 2000;Xing and Zhu, 2000).Summarized from the collected data for model input, the observed average load of point source NH 4 -N into rivers was approximately 4.70×10 4 t year −1 in the Shaying River catchment.The diffuse source contributed 38.57% of the overall NH 4 -N load on average from 2003 to 2005, and this value was slightly higher than the statistical results (29.37 %) given in the official report (Wang, 2011).Moreover, the diffuse source contributions at the stations ranged from 31.72 % (Huaidian) to 47.13 % (Shenqiu).Compared with the diffuse source loads in the individual administrative regions in 2000, the simulated loads tended to increase from 2003 to 2005, except in the Kaifeng region.The yields in the Fuyang and Pingdingshan regions increased at the highest rates.The primary pollution source in the Shaying River catchment was still the point source, but the diffuse source was also an important concern.In terms of spatial variation, the contribution of diffuse source to the pollutant load was high in the upstream and low in the middle and downstream because the point source emission was usually concentrated in the mid-dle and downstream.Therefore, compared with the results in Zhang et al. (2013), the overall simulation performance of NH 4 -N concentration was also improved remarkably by considering the detailed nutrient processes in the soil layers.

Crop yield simulation
The simulated corn yield and its spatial pattern were shown in Fig. 11.The average annual yields were summarized at sub-basin scale and ranged from 0.08 to 326.95 t km −2 year −1 with the average value of 76.84 t km −2 year −1 .The yield of each administrative region was further summarized and compared with the data from statistical yearbooks from 2003 to 2005 (Henan Statistical Yearbook, 2003, 2004and 2005).The high-yield regions were Luohe, Fuyang and Zhoukou in the middle and downstream where the primary land use/cover was the dryland agriculture (93.12, 95.87 and 93.18 %, respectively).The crop yields in the Luohe, Nanyang and Kaifeng regions were well simulated.The total yield was underestimated in the whole basin with a bias of 19.93 %.The discrepancies might be caused by the boundary mismatch between the administrative region and sub-basin, spatial heterogeneities of human agricultural activities and inaccurate cropping pattern used in such huge regions.A high-resolution remote sensing image and field investigation might be helpful to improve the model performance.

Comparison with other models
It is a natural tendency that models grow in complexity in order to capture more interactions of complex water-related processes in the real basins because of more and more available observations and improved accuracies (Beven, 2006).Our proposed model was developed in this direction and tended to benefit integrated river basin management, although the model applicability needs to be further evaluated in different regions.In comparison with most existing models, our proposed model considered all the water-related processes as an integrated system rather than isolated systems for individual processes.
Our model provided competitive simulation results in the Huai River basin (Figs.7-9; Tables 3-5).Several typical models were also applied in this basin, such as SWAT for the monthly runoff and water quality simulation at the regulated stations (Zhang et al., 2013), the SWAT and Xinganjiang models for the daily runoff simulation at the unregulated upstream stations (Shi et al., 2011) and the DTVGM for daily runoff simulation (Ma et al., 2014).Compared with the results of these models, our model generally performed better in the runoff or water quality simulations.In particular, our model performed even better than SWAT at the regulated stations, as more detailed dam regulation rules and soil biochemical processes were considered.For example, the average values of f runoff at the monthly scale decreased from 0.32 (SWAT in Zhang et al., 2013) to 0.15 (our model) at the regulated stations.The average values of f NH4−N decreased from 0.47 (SWAT in Zhang et al., 2013) to 0.27 (our model).Moreover, both the Xinanjiang model and the DTVGM are limited to simulating the flow series at the unregulated or lessregulated stations because they do not consider the dam regulation in their current model frameworks (Shi et al., 2011;Ma et al., 2014).

Equifinality
Until now, our understandings of water-related processes have still been ambiguous, and it is hard to describe all these processes in the real-word systems from strong physical foundations (Beven, 2006;Hrachowitz et al., 2014).Empirical equations are usually adopted to approximate the physical processes with numerous unknown parameters, especially in the large-scale models.A single output variable of models is associated with multiple processes and many parameters.For examples, SWAT contains over 200 parameters (Arnold et al., 1998) and DNDC has nearly 100 parameters (Li et al., 1992).Pohlert et al. (2006) reported that six hydrological and 12 N-cycle sensitive parameters were detected in SWAT-N for the simulation of water flow and N leaching.In the case study, 9 and 14 sensitive parameters of our model were detected for runoff and NH 4 -N simulation, respectively (Table 2).Therefore, due to the large numbers of model parameters and limited observations, most existing models are subject to equifinality, which is more serious if more waterrelated processes are considered or more sub-basins are delineated for the distributed models.
Several strategies would be helpful to alleviate the equifinality, such as field experiments on the physical parameters (Kirchner, 2006), the utilization of more observed processes, multiple evaluation measures for a single predicted component (Her and Chaubey, 2015), parameter regularization and process constraints (Tonkin and Doherty, 2005;Pokhrel et al., 2008;Euser et al., 2013).Moreover, some attempts are made to move away from traditional curve fitting towards more process consistency and efficient model selection techniques (Hrachowitz et al., 2014;Fovet et al., 2015).
For our model, all the independent calibration and validation data sets were specified in Table 1, and most widely used measures of model performances were also provided in the PAT.In the case study, we also employed several observation sources (e.g., runoff and water quality observations at different stations, the diffuse pollution load and crop yield data) and used three measures to evaluate model performance for the individual components (e.g., bias, r and NS).To make full use of the existing data in practice, parameter sensitivity analysis would be an effective way to reduce dimensionality in model calibration and then focus only on the critical processes and parameters that are sensitive to model outputs (van Griensven et al., 2006).Model autocalibration would be efficient to obtain the optimal simulations from numerous samples in multi-dimensional parameter spaces.

Model limitations
It should be noted that our extended model still has several limitations.
1.The mathematical descriptions of groundwater, crop growth processes and agriculture management practices were still inaccurate.The current version focused on the detailed descriptions of hydrological and nutrient cycles in the soil layers and water bodies, and the consideration of dam regulation.Satisfactory performances on water quantity and quality simulation were achieved in our case study.However, the simulations for groundwater, diffuse pollution and crop yield in the agriculture regions could be improved further.The stratification of water impounding in the water quality module should be considered if the high-resolution bathymetric data of dams or lakes are available.
2. High parameterization is an inevitable issue because of its all-inclusive framework.Our model considered the main water-related processes in the hydrological, ecology and water quality subsystems, but numerous processes were still controlled by unmeasurable parameters because of their empirical and/or scale-dependent Hydrol.Earth Syst.Sci., 20, 529-553, 2016 www.hydrol-earth-syst-sci.net/20/529/2016/ nature (Her and Chaubey, 2015).Although the parameter sensitivity analysis and calibration are widely used to handle the high parameterization issue, the equifinality and parameter uncertainty are still inevitable because of the insufficient observations and the complex interactions among different subsystems.

Conclusions
In this study, the TVGM hydrological model was extended primarily to an integrated water system model to address the complex water issues emerging in the basins.The model performance was demonstrated in the Shaying River catchment, China.The model provided a reasonable tool for the effective water governance by simultaneously simulating several indicative components of water-related processes including the hydrological components (e.g., runoff, soil moisture, evaporation and plant transpiration, water storage in the dams and sluices), water quality components (e.g., diffuse pollution source load, water quality concentrations in water bodies) and ecological components (e.g., crop yield), which could be calibrated if observations were available.The case study showed that the simulated runoffs at most stations fitted the observations well in the highly regulated Shaying River catchment.All the evaluation criteria were acceptable for both the daily and monthly simulations at most stations.This model simulated the discontinuous daily NH 4 -N concentration well and properly captured the spatial patterns of diffuse pollution load and corn yield.
Owing to the heterogeneity of spatial data in large basins and insufficient observations of individual subsystems, not all the results were acceptable and several processes were still not well calibrated (such as low flow events, diffuse pollution source load and crop yield).More available data and improved data quality will reduce the model uncertainty and equifinality problem, especially the higher-resolution data for surface conditions, water quality, agricultural management and socio-economic data.The model would be improved by further considering more accurate human activities in the agricultural management, calibrating multiple components by multi-objective optimization and model uncertainty analysis because of the interactions and tradeoffs among different processes.The over-parameterization and the reasonable prior parameter conditions should also be treated carefully in applications.Advanced analysis technologies would benefit the future model development, such as model selection techniques, parameter regularization.Moreover, an easily used operational software package can broaden the model's applications in different regions.More case studies are needed to further demonstrate its applicability.and k 2 are the specific decomposition rates of labile faction and resistant fraction, respectively (day −1 ).
The NH 4 amount (FIX NH4 , kg ha −1 ) absorbed by clay and organic matter is estimated by FIX NH 4 = 0.41 − 0.47 • log(NH 4 ) • (CLAY/CLAY max ) , (B4) where NH 4 is the NH + 4 concentration in the soil liquid (g kg −1 ).CLAY and CLAY max are the clay content and the maximum clay content, respectively.The Supplement related to this article is available online at doi:10.5194/hess-20-529-2016-supplement.

Y
. Y. Zhang et al.: Integrated water system simulation velopment of water-related sciences, computer science, Earth observation technologies and the availability of open data.

Figure 2 .
Figure 2. The flowchart of HCM and the interactions with other modules.

Figure 3 .
Figure 3.The flowchart of SBM (a) and CGM (b) in the ecological part and the interactions with other modules.

Figure 4 .
Figure 4.The flowchart of SEM (a), OQM (b) and WQM (c) in the water quality part and the interactions with other modules.
BOD and COD.Point pollution sources are also consid-Hydrol.Earth Syst.Sci., 20, 529-553, 2016 www.hydrol-earth-syst-sci.net/20/529/2016/ ered.Point sources are directly added to the surface water in the model according to their geographic positions.Common point sources are urban water treatment plants and industrial plants.

Figure 5 .
Figure 5.The flowchart of PAT and its interactions with other modules.

Figure 6 .
Figure 6.The location of the study area (a) and the digital delineation of the sub-basin, point source pollutant outlets, rural population (b), animal stock (c) and fertilization (d).

Figure 7 .
Figure 7.The daily runoff simulation at all stations.

Figure 8 .
Figure 8.The cumulative distributions of simulated and observed daily runoff at all stations.

Figure 11 .
Figure 11.The spatial pattern of corn yield at the sub-basin and regional scale in the Shaying River catchment.

Table 1 .
The data sets and their categories used in the model.

Table 2 .
Sensitive parameters, their value ranges and relative importance for runoff and NH 4 -N simulations.

Table 3 .
Runoff simulation results for regulated and less-regulated stations.