Integrated hydrological modeling of the North China Plain and implications for sustainable water management

Groundwater overdraft has caused fast water level decline in the North China Plain (NCP) since the 1980s. Although many hydrological models have been developed for the NCP in the past few decades, most of them deal only with the groundwater component or only at local scales. In the present study, a coupled surface water–groundwater model using the MIKE SHE code has been developed for the entire alluvial plain of the NCP. All the major processes in the land phase of the hydrological cycle are considered in the integrated modeling approach. The most important parameters of the model are first identified by a sensitivity analysis process and then calibrated for the period 2000–2005. The calibrated model is validated for the period 2006–2008 against daily observations of groundwater heads. The simulation results compare well with the observations where acceptable values of root mean square error (RMSE) (most values lie below 4 m) and correlation coefficient ( R) (0.36–0.97) are obtained. The simulated evapotranspiration (ET) is then compared with the remote sensing (RS)-based ET data to further validate the model simulation. The comparison result with a R2 value of 0.93 between the monthly averaged values of simulated actual evapotranspiration (AET) and RS AET for the entire NCP shows a good performance of the model. The water balance results indicate that more than 70 % of water leaving the flow system is attributed to the ET component, of which about 0.25 % is taken from the saturated zone (SZ); about 29 % comes from pumping, including irrigation pumping and non-irrigation pumping (net pumping). Sustainable water management analysis of the NCP is conducted using the simulation results obtained from the integrated model. An effective approach to improve water use efficiency in the NCP is by reducing the actual ET, e.g. by introducing watersaving technologies and changes in cropping.


Introduction
The use of groundwater as a resource for water supply to urban areas and for irrigation has increased dramatically during the past decades, while in many areas of the world there is still a large unused potential for increased groundwater exploitation (Zektser and Everett, 2004).In some areas of the world groundwater levels of major aquifers are rapidly declining as an effect of unsustainable pumping and insufficient recharge.This is for instance the case in West Java, Indonesia; Dhaka, Bangladesh; and Beijing, China (Braadbaart and Braadbaart, 1997;Hoque et al., 2007;Zhou et al., 2012).Large-scale hydrological models have become increasingly popular in assisting water resources management of large aquifer systems (Casper and Vohland, 2008;Goderniaux et al., 2009;Branger et al., 2010).However, great challenges exist in developing and validating such models mainly due to data issues such as unevenly distributed monitoring sites or simply insufficient data (van der Linden and Woo, 2003;Cao et al., 2006;Ferguson and Maxwell, 2010).As a result, the reliability of the model simulations is in question.

H. Qin et al.: Integrated hydrological modeling of the North China Plain
Remote sensing (RS) data by satellite may have the potential to provide an ideal solution to the problem of data scarcity for large-scale hydrological models.RS data are normally spatially distributed and with large coverage.This feature is especially beneficial in areas with sparse observational data.However, there is a lack of experience on implementing RS data in large-scale distributed hydrological models since the method of effectively using such data to constrain or validate hydrological models is rather difficult (Chen et al., 2005;McMichael et al., 2006;Stisen et al., 2008).Therefore, the question remains how we can obtain reliable model simulations by taking advantage of the added value from the extra information RS can provide.
The North China Plain (NCP) is one of the most densely populated areas in the world.Due to the rapid economic growth in the past few decades, the over-exploitation of groundwater resources has constantly raised the alert.Numerous studies, most of which solely focus on the groundwater component and employ MODFLOW as their modeling tool, have been completed to analyze the NCP's hydrological system and groundwater resources utilization, such as Jia et al. (2002), Foster et al. (2004), Mao et al. (2005), Xu et al. (2005), Nakayama et al. (2006), Wang (2006), Xu (2006), Liu et al. (2008), Wang et al. (2008), Shao et al. (2009), Zhang et al. (2009), Lu et al. (2011), Zhou et al. (2012) and Cao et al. (2013).Groundwater recharge in these models is oversimplified, estimated as a fraction of the precipitation without considering the dynamics of the evaporation or the existence of irrigation.Kendy et al. (2003) pointed out that the recharge fraction depends on the quantity and temporal distribution of precipitation and irrigation.The oversimplified treatment of recharge introduces uncertainty in evaluating the groundwater recharge and hence in closing the water balance.The Institute of Remote Sensing Applications (IRSA), Chinese Academy of Sciences provided estimates of ET over the NCP using a combination of data from meteorological stations and RS data (Wu et al., 2012).In order to make full use of such ET data and hence make a more reliable assessment of the water balance for the NCP aquifer, it is required to use an integrated hydrological model with a coupled groundwater-surface water approach (Henriksen et al., 2003;Doummar et al., 2012;Li et al., 2012;Stisen et al., 2012).In the present study, the recharge component is treated as a seasonal variable and is distributed based on the formulations of evapotranspiration and unsaturated flow.Precipitation, irrigation, potential evapotranspiration, and crop and soil properties are also considered in the model to calculate the recharge.
The objectives of the present study are to develop an integrated hydrological model for the NCP aquifer system using the MIKE SHE modeling tool and to analyze and understand the water balance components of the NCP.The model is used to assess the potential of using the RS-based evapotranspiration in order to improve the water balance calculations and to discuss the perspectives of using an integrated hydrological model to support sustainable water resources management in the NCP.

The North China Plain
Geographically, the North China Plain (NCP) refers to the alluvial plain north of the Yellow River, south of the Yan Mountain and east of the Taihang Mountain.Figure 1 shows the region of the NCP, spanning from 32 • N to 40 • N and from 114 • E to 121 • E. It includes all districts of Beijing and Tianjin Municipality, most of Hebei Province and the northern area of the Yellow River Basin of Shandong and Henan Province.The total area of the NCP is 140 000 km 2 and it is one of the most densely populated regions in the world, with an average population density of approximately 800 persons km −2 .In the year 2000, the total population of the NCP was around 130 million.It is the center of China's political, economic and cultural development.Due to its subhumid to semi-arid climate condition, the NCP is very prone to drought, which stresses sustainable water resources management of that region (Liu et al., 2008;Zheng et al., 2010;Qin et al., 2012).
The NCP aquifer is divided into four main hydrogeological zones from the Taihang Mountains on the west to the Bohai Sea on the east (Fig. 1); these zones are based on the geomorphology of paleochannels.According to their sedimentary character and relative geographic positions, the region consists of a piedmont plain, an alluvial fan, an alluvial plain (central plain), and a coastal plain (Wu et al., 1996;Cao et al., 2013).
The climate in the NCP is continental semi-arid, with average annual temperature of 12-13 • C. The mean annual precipitation of 1951-1995 was 554 mm, while the seasonal distribution of precipitation is uneven, with about 75 % of its precipitation occurring throughout the summer flood season from July to August (Qin et al., 2012).This climate type is the reason for frequent flooding and drought events in the NCP region.There have been 135 flooding years and 140 drought years in the NCP in the last 500 yr (Huai River Commission, 1998).The water resource per capita in the NCP is 501 m 3 , which is only 23 % of China's average (Xia et al., 2007).Groundwater has become an indispensable factor to support the economic and social sustainable development of the NCP.Currently, the groundwater resources have accounted for more than 70 % of the total water supply in Beijing, Shijiazhuang and many other cities in the NCP.Overexploitation of the groundwater resources is very common, especially during the past few decades, where the acceleration of industrialization in China took place.Moreover, the industrialization and growth of the population are unsustainable and the economic development in the area is likely to be limited due to the scarcity of water in near future (Qin et 42 Location and topography of the North China Plain, spanning from 32°N to 40°N and 114°E to 121°E (Cao et al., 2013).( 1 (Cao et al., 2013).
(1) Locations of meteorological stations and observation wells: the wells marked with rectangles are the selected wells used to show the calibration and validation results of groundwater heads; (2) boundaries of geomorphologic zones: (I) piedmont plain, (II) alluvial fan, (III) alluvial plain, and (IV) coastal plain (Wu et al., 1996;Cao et al., 2013); (3) locations of capital cities (Beijing, Tianjin, Shijiazhuang and Jinan) and some other cities that are included and considered in the NCP model.al., 2012).Agriculture accounts for 85 % of the total water consumption in the NCP (Xu et al., 2005).Winter wheat and summer maize are the main crops in the NCP rotation system, accounting for about 80 % of total cultivated area and 90 % of grain yield in the NCP (National Bureau of Statistics of China, 2010;Shu, 2010).
Table 1 lists the types and sources of the data used in the NCP model setup.Data defining the driving forces for the water flow is described in the following sections.

Climate data
In the NCP model, climate data include precipitation rate and reference evapotranspiration.The spatial distribution type of precipitation rate is fully distributed, whereas that of the reference evapotranspiration is station-based.The precipitation data applied in the model are the 3B42 daily product from the Tropical Rainfall Measuring Mission (TRMM) of NASA Science (http://trmm.gsfc.nasa.gov)with a resolution of 0.25 • .The TRMM 3B42 product is a merged product containing satellite information from both microwave and IRsensors and is calibrated against rain gauge measurements on the monthly scale (Huffman et al., 2007).The precipitation data are resampled to 10 km resolution and re-projected to match the MIKE SHE model grid using bilinear resampling (using GDAL-Warp).The data of reference evapotranspiration are generated using data from 19 meteorological stations (locations are shown in Fig. 1; 10 in Hebei Province, 5 in Shandong Province, 2 in Henan Province, 1 in Beijing Municipality and 1 in Tianjin Municipality) provided by the Chinese Meteorological Data Service (CMDS) and calculated using the Penman-Monteith equation (Monteith, 1965).The station-based reference evapotranspiration data are distributed using Thiessen Polygons.Table 2 lists the information of the meteorological stations, including station name, elevation, latitude, longitude, average yearly reference evapotranspiration, mean temperature and which municipality or province this station belongs to.

Land use data
The land use map (listed as Fig. 2a) is created from MODIS NDVI products from 2003.The classification is made from the seasonal variation in NDVI using the method proposed by Kang et al. (2007) and Shu et al. (2012).It is assumed that the generated land cover classification can represent the conditions throughout the whole simulation period.The land cover classification is divided into nine classes: water, resident, wheat and corn_piedmont, rock, fruit tree, cotton_central, cotton_coastal, wheat and corn_central, and wheat and corn_coastal, where "piedmont" refers to piedmont district, and "central" and "coastal" refer to the central and coastal plains, respectively.The land use type wheat and corn is divided into three types, while the type cotton is divided into two types.One reason for doing so is that the vegetation parameters (e.g.leaf area index, crop coefficient and root depth) are different in piedmont plain, alluvial plain and coastal plain (Ko et al., 2009;Piccinni et al., 2009).Another reason is that even the same crop needs different amounts of irrigation water for different plains in the NCP (Ma et al., 2011).The vegetation parameters are assessed on the basis of experience or collected from literature.In this study, the leaf area index (LAI) and root depth take the same values for a certain vegetation type, while crop coefficients (K c ) and irrigation water requirements have different values for the same vegetation type.These are required input data for the model setups.Table 3 lists the typical value ranges of LAI, root depth and K c for each crop in the NCP.
The typical crop rotation in the NCP is winter wheat and summer corn.The growth of winter wheat is described by six stages from October to early June (approximately 210 days): sowing, wintering, revival, jointing, flowering, and grain-filling.And from June to September the summer corn is grown as described by four stages: sowing, emergence, jointing-flowering, and grain-filling (Shu, 2010).The growing season of cotton and fruit trees is the same, which is from May to October.Irrigation is applied according to the growing stages of different crops.The irrigation frequency and yearly irrigation water amount for these crops are listed in Table 3.

Pumping data
Groundwater has been used for industrial and domestic, and irrigation purposes in the NCP.The pumping for domestic and industrial water use and the pumping for irrigation are estimated separately.The former is approximated using a system dynamics (SD) model incorporating multiple sources of data (Qin et al., 2012); while the latter is simulated using the irrigation module embedded in the MIKE SHE code.
In the NCP, farmers use shallow groundwater wells for irrigation, where the locations of the wells may not be mapped exactly.No information exists on the number of pumping wells used for irrigation.In a large regional-scale model like the NCP, each grid cell could thus contain many shallow groundwater wells.Therefore, the "Shallow Well source" in MIKE SHE has been selected to simply extract water for irrigation from the same cell where it is used, without having to know the exact coordinates of the wells.By specifying this option, one well is placed in each cell of the irrigated area.The water application method for irrigation is sprinkler, which means that irrigation in the model is added to the precipitation and the amount and timing of irrigated return flow is explicitly calculated by the model.Irrigation water amount is estimated according to different crop types.Different irrigation frequency and water amount are applied to different crops.Take winter wheat for instance; there are five irrigation times in one year and each time winter wheat is irrigated with a certain amount of water.The yearly irrigation water of winter wheat can be computed by multiplying the irrigation frequency and irrigation water needed each time.Then the total irrigation water is the sum of water amounts of all the crops in this area.
Most of the water supply to the cities in the NCP is abstracted from groundwater including both the shallow and deep aquifer systems.We consider groundwater abstraction for domestic and industrial use in all the 21 cities (Beijing, Tianjin, Qinhuangdao, Tangshan, Shijiazhuang, Handan, Xingtai, Baoding, Langfang, Cangzhou, Hengshui, Jinan, Dongying, Binzhou, Dezhou, Liaocheng, Anyang, Hebi, Xinxiang, Jiaozuo and Puyang) in the entire NCP aquifer (Qin et al., 2012).The capital cities (Beijing, Tianjin, Shijiazhuang and Jinan), which are the most important industrial and domestic areas, and some other cities are shown in Fig. 1.Abstraction data of these cities are provided by Qin et al. (2012) and the pumping wells are assumed to pump water from the deep aquifer zone (layer 2 and 3, or aquifer units III and IV).
There are also so many wells pumping groundwater for domestic and industrial purposes that they cannot be mapped exactly.To represent the pumping activities in the NCP, a random-point approach and the concept of "virtual" well have been adopted in this study.The random-point approach is applied to the 21 aforementioned cities in the NCP.To consider the different recycling ratios of urban and rural water uses, the domestic water demand is divided into urban domestic and rural domestic (including rural population and  livestocks) water demand and a fraction F is introduced here by simply reducing the groundwater pumping amount to reflect the recycling effect, as the pumped water here disappears from the system.Specifically, the groundwater pumping for non-irrigation use is formulated based on the hypotheses from Qin et al. ( 2012): (1) collect time series data of total urban domestic water demand UDWT (i, t), total rural domestic water demand RDWT (i, t) and total industrial water demand IWT (i, t) of the ith city, respectively, where 1 ≤ i ≤ 21, and t is the simulation time.The total nonirrigation water demands are estimated based on the population and industrial output using the SD method (Qin et al., 2012).( 2) Randomly generate N i points in the ith city area, respectively, where 1 ≤ i ≤ 21 and 10 ≤ N i ≤ 25.The random point's number of each city depends on the ranking of relative importance of that city.For example, there are 25 points in Beijing City and 10 points in Jiaozuo City.Each random point is assumed as a "virtual" pumping well in that area and the same pumping rate is applied on each "virtual" well.(3) Collect the consumption rates of urban domestic water demand (F UD ), rural domestic water demand (F RD ) and industrial water demand (F I ) from literature.The consumed water is the one that disappears from the hydrological system.The surface runoff into the sea is another source of water that disappears from the system.The average runoff into the sea is about 10.6 mm yr −1 (Yu, 2008).This amount of water is assigned to the urban domestic sector and industrial sector equally.That is to say, both the urban domestic sector and industrial sector will "consume" an extra 5.3 mm yr −1 of water.(4) For the ith city, calculate the pumping rate for each well using Eq. ( 1):  Table 3.Typical values for LAI, K c and root depth for different crops (Liu et al., 2002;Kendy et al., 2004;Zhang et al., 2004;Shu et al., 2012) and irrigation frequency and total irrigation water amount for these crops.where UDW(i, t), RDW(i, t) and IW(i, t) are the urban domestic, rural domestic and industrial water demand pumping rate for each well in the ith city, respectively; α UD and α I are the fractions of 5.3 mm yr −1 to the total urban domestic and industrial water demand, respectively.A total number of 870 "virtual" wells (including urban domestic, rural domestic and industrial pumping wells) are generated and distributed into the NCP aquifer after these four steps.

ETWatch data
ETWatch has originally been developed as an operational software system for regional ET monitoring, as per request   , 2005, andCao et al., 2013).Aquifer units I and II make up layer 1 and also shallow aquifers in the model, while aquifer units III and IV make up layer 2 and 3, respectively; aquifer III and IV are deep aquifers here.
of the Hai Basin commission (Wu et al., 2008a;Xiong et al., 2010a, b).ETWatch is an integration of remote sensingdriven surface energy balance models and a Penman-Monteith (P-M) module used during cloudy conditions (Wu et al., 2011a(Wu et al., , b, 2012)).It is an operational processing chain starting from data pre-processing to application products, utilizing spatial information of climate, soil type, land use, vegetation cover, digital elevation and remote sensed land surface parameters.The surface energy balance algorithm for land (SEBAL) and the surface energy system (SEBS) are integrated into ETWatch to estimate surface fluxes during clear-sky conditions (Wu et al., 2008b).The ETWatch data are provided by IRSA of the Chinese Academy of Sciences and used in the ET comparison section of our study to validate the model simulation (Wu et al., 2012).
3 Model setup and calibration

Modeling tool
The modeling tool used in this study, MIKE SHE, is developed by DHI (http://www.dhigroup.com).It was originally developed as the SHE (Système Hydrologique Européen) model (Abbott et al., 1986a, b) along with the British Institute of Hydrology and SOGREAH (France).MIKE SHE is a deterministic, distributed, physically based hydrologic modeling system that covers the major processes in the hydrologic cycle and includes process models for evapotranspiration, overland flow, channel flow, unsaturated flow, groundwater flow, and their interactions.A finite difference approach is used to solve the partial differential equations describing the processes of overland (two-dimensional Saint-Venant equation), channel (Kinematic Routing solver), unsaturated (2-layer water balance method) and saturated subsurface flows (three dimensional Boussinesq equation), whereas analytical solutions are used for describing interception and evapotranspiration (Refsgaard and Storm, 1995;DHI, 2012a, b).Since the precipitation is small in winter and the data are not available for the NCP model, snowmelt is not activated in the model simulation.Other major components included to describe the hydrological system are the irrigation and the groundwater pumping.

Model discretization
The integrated MIKE SHE model of the NCP covers the period from 1 January 2000 to 31 August 2008, a total of nearly nine years.Due to the different timescales in different flow systems, the time step of groundwater flow is one day while that of the river flow is five minutes.The grid covering the modeling area is discretized to 300 columns and 350 rows with a uniform cell size of 2 km × 2 km, a total simulation area of 140 000 km 2 .Based on the first hydrogeological survey conducted in the 1970s by the Hebei Bureau of Geology and Minerals, the Quaternary aquifers in the NCP were divided into four major aquifer units (Chen et al., 2005;Cao et al., 2013), which are shown in Fig. 3 as units I-IV.The geologic units corresponding to these aquifers are Holocene series (Qh 4 ) and late (Qp 3 ), middle (Qp 2 ) and early (Qp 1 ) Pleistocene series (Cao et al., 2013;Foster et al., 2004;Zhang et al., 2009).The model grid is discretized into three layers for the aquifer system (Cao et al., 2013).The first layer of the grid includes the first two aquifer units (I and II in Fig. 3), while the second and third layer correspond to aquifer unit III and IV.The Shuttle Radar Topography Mission (SRTM)  ( Rabus et al., 2003) elevation data at 90 m resolution have been aggregated to the model grid and used as the top elevation (Cao et al., 2013).

Rivers and lakes
There are many rivers in the entire NCP area and the five major rivers -Luan River, Yongding River, Hutuo River, Zhang River and Wei River -are considered in the NCP model (Fig. 4).These rivers are selected due to data availability and the fact that the unconsidered rivers in the NCP have run dry after the year 2000, which is the start year of the simulation period.The network, cross section and boundary of the rivers are input data of this section.These data are provided by the China Institute of Water Resources and Hydropower Research (IWHR).The geometry of each river branch is specified in terms of cross sections.Cross sections and data are important to both conveyance capacity and storage capacity at different reaches of river system.Two types of river boundaries are considered for each river branch in the study: inflow at the upstream and water level at the downstream.The latter boundary is a constant one due to the data availability of the rivers.Each river exchanges water with groundwater zone through the coupling reaches.Each coupling reach is interpolated to a MIKE SHE river link, which is located at the nearest boundary between grid cells (DHI, 2012a, c).

Unsaturated zone
Soil types and their related properties are very important for the unsaturated flow.Figure 2b is the soil type map of the NCP, with loam as the dominant type and silty loam and sand as secondary types.Sandy soil is mainly located along Hutuo River and in other isolated patches within the model area.Other sediment types, like silt loam is only found at few locations (Shu, 2010).The 2-layer water balance model for the unsaturated zone requires the following data of each soil type (Shu et al., 2012): water content at saturation (θ sat ), water content at field capacity (θ fc ), water content at wilting point (θ wp ) and saturated hydraulic conductivity (K s ).The values of the three soil type parameters are listed in Table 4.Additional parameters are required for MIKE SHE in the 2layer water balance method, which is based on a formulation presented in Yan and Smith (1994), to calculate the actual ET.The ET surface depth, which equals the thickness of the capillary zone, is required.This thickness is added to the root depth to define the thickness of Layer 1 in the 2-layer water balance ET method.In unsaturated soils, capillary action can lead to saturated conditions existing some distance above the water table.ET will not decrease the saturation, but draw water directly from the water table due to capillary action.If the water table is at or above the bottom of the root zone then ET will be removed from saturated zone (SZ) at the full rate.
As the water table falls below the root zone, ET from SZ will be reduced linearly until the water table is at a distance below the root zone equal to the ET surface depth -at which point the ET from SZ becomes zero (DHI, 2012a).The value of the ET surface depth is set to 0.2 m in the NCP model to represent the loamy and sandy soil types.

Saturated zone
The groundwater zone is conceptualized into 3 layers describing a multi-layer, heterogeneous and anisotropic 3-D hydrogeological model of the NCP aquifers.The layers extend to a depth of more than 500 m below the land surface.
The hydrogeological conceptualization proposed by Cao et al. (2013) 2c and d) for the reasons that it represents the geological and hydrogeological characteristics of the NCP aquifer.The geological units (and their parameters) adopted from Cao et al. (2013) represent categorized zonal index using a continuous K field, therefore they do not have hydrogeological meanings.The same parameters (horizontal conductivity K xx and K yy , vertical conductivity K zz , specific yield S y , and storage coefficient S s ) are used in the same unit.The hydraulic conductivity and storage coefficients were estimated using pumping test data and their values were summarized for the different lithologies across the NCP by Chen (1999) and Zhang et al. (2009).Hydraulic conductivities, specific yields and storage coefficients were mapped across the NCP by the China Geological Survey (CGS) (Zhang et al., 2009).Layer 1 is a shallow aquifer zone while layers 2 and 3 are deep aquifer zones.As suggested by Cao et al. (2013), layers 2 and 3 are described by the same distribution of hydrogeological units (Fig. 2d).Initial potential heads of 1 January 2000 for layers 1, 2 and 3, which are used as the initial conditions in the model simulation, are adopted from Cao et al. (2013).Three types of boundary conditions are assigned to layer 1: fixed head to zero, fixed head to initial value, and zero-flux.The first one is used at the northeastern area of the model, where the land intercepts with the Bohai Gulf, to represent the average sea level.The second one is used at the south and southeastern parts of the model where it is bounded by the Yellow River.The last one is used to the western boundary, which is the interface between the Taihang Mountain and the piedmont plain.For layer 2, two boundary conditions are used: the northeastern part, which is the interface with the Bohai Gulf, is assigned a zero-head condition; the rest of the area is assigned a zero-flux condition, which is also applied for the whole boundary for the third layer.

Model calibration
The primary aim of the calibration process is to obtain a set of model parameters that provide a satisfactory agreement between simulation results and field observations.The NCP model is calibrated using AutoCAL as an optimization tool, which is included in the MIKE SHE software package (Madsen, 2003).AutoCAL applies the shuffled complex evolution (SCE) optimization algorithm (DHI, 2012d), which is a global optimization algorithm that combines various search strategies, including (i) competitive evolution, (ii) controlled random search, (iii) the simplex method, and (iv) complex shuffling (Madsen, 2003).
The calibration targets are groundwater heads, which are very important in the model simulation and will influence saturated zone flow and actual ET.There are 226 observation wells (locations are shown in Fig. 1) distributed in the study area that are used in the calibration procedure.Calibration results of 9 selected observation wells (shown as selected wells with names in Fig. 1) are shown in Fig. 5.The criteria for the selection of these nine observation wells are that they need to have complete time series for the entire model simulation period and that they are less disturbed by the local pumping activities.Table 5 lists the information of these selected observation wells.The information presented here includes the geographical locations, the ground depth and the screen depth where the groundwater head is measured.The simulation period is divided into two stages: the calibration stage from 2000 to 2005 and the validation stage from 2006 to 2008.
The root mean square error (RMSE), which is an aggregated measure that includes both the bias and the dynamical correspondence (Anderson and Woessner, 1992) and can be calculated by Eq. ( 2), is the objective function minimized in the auto-calibration process of the NCP model.Furthermore, there are another four quantitative performance criteria as defined in Eqs. ( 3) to (6): ME (mean error), MAE (mean absolute error), STD res (standard deviation of residuals) and R (coefficient of correlation).
where O i is the ith (i is the daily index) observed value, S i is the ith simulated value, O is the average observed value, S is the average simulated value, N is the total number of corresponding measurements and simulations.

Calibration and validation results
Prior to the calibration, a sensitivity analysis is performed using AutoCAL to identify the most sensitive model parameters to be fine-tuned in the succeeding parameter optimization.The local sensitivity analysis method (DHI, 2012d), which provides the sensitivity of the model parameters around a specified parameter set and hence gives information about the importance of the parameters only at that location in parameter space, is adopted in our model.Equation ( 2) is used as the objective function in the auto-calibration process and equal weight is set for the parameters.
The parameters related to the unsaturated zone are adopted from Shu (2010) and are therefore not included in the calibration procedure.The parameters selected as the most sensitive parameters and then used in the calibration stage are listed in Table 6.The final values of the parameters, i.e. the optimal parameters after the calibration, are also listed in Table 6.The most sensitive parameters are specific yield for geological units and horizontal conductivities, where the latter are log-transformed.Some parameters are tied to each other in order to keep the number of free parameters as low as possible (Table 6).The calibration performance criteria calculated of each observation well for the calibration and validation periods are listed in Table 7.At the calibration stage, the performance criteria suggest reasonable calibration results for the model.For most observation wells, the correlation coefficient is above 0.75, which indicates good agreement between the observed and the simulated groundwater heads.The values of RMSE of most wells are between 0.86 m and 3.51 m, which are considered small for such a large simulation area.The deficiencies of RMSE, ME and MAE are likely to be caused by the aggregation and simplification of parameters for the groundwater flow system.At the validation stage, except for wells #116 and #140, the correlation coefficient for most wells are very high, which shows that the calibration is rather successful.The performance result for the total 226 wells (with one average value on behalf of that well) is also listed in Table 7.The ME, MAE and RMSE values are 1.43 m, 13.11 m and 19.44 m, respectively.These values are considered small relative to the maximum groundwater level variations in the primary cones of depression of about 60 m and 70 m in the shallow and deep aquifer zones, respectively (Cao et al., 2013;Fei et al., 2009).Overall the results of the Table 5. Information of the selected groundwater observation wells used in the calibration and validation processes (Cao et al., 2013).The information includes the location information, ground elevation and monitor depth information that determines how deep the observed heads are measured in. Figure 5 lists the comparison of simulated and observed groundwater heads of selected observation wells distributed in the NCP for the calibration period and validation period.The simulated groundwater dynamics fits well with the observations during both periods.A declining trend is observed in all the simulated and observation wells.For a number of wells, the groundwater head dynamics is not reproduced by the simulation.This phenomenon may be explained by the influence of local pumping activities.The intensive pumping activities take place at the same model grid, which may have a smoothening effect on the shape of the simulated head curve; or there may be lack of knowledge on locations of all wells.There are discrepancies in the spatial supporting scales between the simulated and observed groundwater heads.The associated scale uncertainty arises mainly for two reasons: first, the simulated groundwater head usually represent the values at the center of the model grids; however it is quite seldom that the boreholes are located precisely at those places.Second, the model is not able to account for the local heterogeneity at the scale smaller than the model grid.As a result, both the model simulation and the field data could be reasonably correct at their own supporting scales, but the calculated uncertainty (e.g.RMSE) is not minimized to zero.Furthermore, the method used to identify the locations of the groundwater abstraction wells in the cities, as described in Sect.2.4, could also introduce considerable uncertainty.In wells #66 and #116, the simulated curve and the observed curve are in close agreement, which are reflected by the high correlation coefficients and the model efficiency values.The validation results in Fig. 5 further support the reliability of the calibrated hydrological model of the NCP. Figure 6 lists the comparison of averaged groundwater head between the observed and simulated values for the total 226 wells.Each point in Fig. 6 represents one well with average simulated/observed values.This comparison illustrates the level of general reliability of the model over the entire NCP area.

Comparison of actual evapotranspiration
For comparison with the model results, estimates of actual evapotranspiration (AET) are obtained from the remote sensing-based ETWatch data provided by IRSA (Wu et al., 2011a(Wu et al., , b, 2012)).The ETWatch data are available for the period January 2002 to August 2008.
Figure 7 shows a comparison of the monthly AET averaged over entire model domain, excluding an area to the southeast where the ETWatch data are not available.Our model offers a good description of the seasonal dynamics of AET, as reflected in Fig. 7, including the double peak for each year.The double peak is determined by the winter wheat and summer corn rotation type of the NCP agriculture.The model-simulated AET fits very well with the RS based AET, with a R 2 value equal to 0.93, RMSE value equal to 10.7 mm month −1 and MAE value equal to 8.5 mm month −1 .The model-simulated ET values tend to be underestimated compared to the ETWatch estimates especially for higher values.
However, the simulated mean monthly AET are not always in good agreement with the ETWatch values when considering the spatial distribution across the NCP, especially during the period with high AET.The discrepancy is not surprising due to the very different parameterizations of ET in the two methods.Figure 8   features can be found in both estimates of AET, but there are also significant differences.For the model-simulated values, the northern part of the model domain lacks some information on the land cover, causing the model to overestimate the AET rates for urban areas such as the Beijing area.Another notable difference is the occurrence of an area with very high AET values along the coastline in the ETWatch product.This is not the case for the model simulations where this area is dominated by low AET values.
The spatial patterns of AET stratified by the land cover information used in the model are analyzed (described in Sect.2.3).Eight classes (the class type Rock is ignored because it has no data and no effect on the model) are considered to represent different land covers (e.g.agricultural areas and urban areas), as well as different irrigation schemes applied in different regions (e.g.there are several cotton classes; see Fig. 2a).Overall, most classes performed similarly to what is shown in Fig. 7, with linear regressions of R 2 values > 0.92 and slopes around 0.85.The main difference between the simulated AET values and the ETWatch estimates is that they show different spatial patterns in each land use class.For ETWatch, there are much larger differences between the different land cover classes than the case for the model-simulated values.As no independent validation data are available, it is not possible to further evaluate the results of the two ET products, but it is clear that the modelsimulated ET has more homogenous fields compared to that of the ETWatch.This is partly due to the differences in parameterization of the land cover in the two methods.Figure 9 shows the mean annual AET for the two different methods.We have also added the variance in terms of the standard   deviation, calculated based on the average values for each of the eight classes considered.For all years except 2006, the standard deviation is higher for the ETWatch data than for the model estimates.
The two approaches in estimating AET are fundamentally different.However, the AET comparison indicates that to some extent the vegetation parameters of our model are reasonable since the magnitude and timing between the two AET estimates are comparable based on the monthly data.As AET is a very important component in NCP's water balance analysis, this comparison of AET is a useful test of the NCP model.the entire NCP.The discharges of the rivers decrease during the simulation period, which is reasonable since most of the rivers in the NCP have run dry in recent years.The discharges are varying within the year, with a peak during the summer flood season.Due to the lack of time series observations, quantitative comparison between simulated and observed stream flows is not presented.However it could bring new insights to the model calibration using stream flow data, when the observed stream flow data become available in the future.

Water balance analysis
Water balance is an important and essential result of the model simulation, which gives us helpful information on available water resources availability and use in the study  area.Water balance analysis is a useful tool in the process of the NCP's sustainable water management.Figure 11 shows the yearly averaged water balance of the entire NCP aquifer system for the whole simulation period while Table 8 shows the annual water balance components for the simulation period.The main components included in the total water balance of the NCP model are precipitation, irrigation, evapotranspiration, and pumping, where the first two are water inflow of the groundwater zone and the latter two are water outflow of the groundwater zone.In general, the total inflow is 737 mm yr −1 while the total outflow is 801 mm yr −1 .This deficit is balanced by a storage depletion of 8 mm yr −1 for the unsaturated zone and 55 mm yr −1 for the saturated zone (Fig. 11).Wang (2006) developed a three-dimensional numerical groundwater model using Groundwater Modeling System (GMS) for the NCP and the calculated groundwater depletion is 52 mm for the years 2002 and 2003.Wang et al. (2008) evaluated NCP's groundwater depletion with a groundwater depletion of 51 mm from January 2002 to December 2003.Compared to these results, the groundwater depletion during this period of our model is 67 mm, which is comparable with Wang (2006) and Wang et al.'s (2008) results.Xu (2006) used GMS to simulate the groundwater flow of Beijing Plain (part of the NCP) and the water deficit for the year 2000 is 42 mm.However, the water deficit simulated by this study for the same period is 53 mm.All these comparison demonstrates the reliability of our model in the hydrological system modeling of the NCP from the aspect of water balance.
Most of the water (70 %) entering the model area leaves the hydrological system through evapotranspiration, of which 99.75 % comes from the unsaturated zone and 0.25 % comes from the saturated zone, while 29 % leaves the system through pumping.This shows an obvious deficit in the total water balance and Fig. 11 shows that the total water going out of the system is more than that entering the system.Except for some wet years, such as 2003, 2004and 2008, Table 8 , Table 8 shows that the actual ET is higher than the precipitation, which indicates an unsustainable use of the water resources in the NCP.Furthermore, the total evapotranspiration is larger than the total precipitation for the entire simulation period, a corresponding result of which is the decline of saturated zone and unsaturated zone storage.The total decline in the saturated zone and unsaturated zone storage during the entire simulation period are 480 mm and 65 mm, respectively.A declining trend of the storage indicates an unsustainable groundwater exploration in the study area (Fig. 12).The groundwater storage change follows an annual cycle with a decrease during the dry winter and increase during the wet summer.As noticed in Table 8, the values of UZ storage change in 2003, 2007and 2008and SZ storage change in 2003 are positive, which indicates an increase in storage due to the high precipitation in these years.The high annual precipitation in 2003 causes an overall increased groundwater storage which is also reflected in an associated rise of the groundwater table (Fig. 5).Compared with other components, the base flow from and to the rivers is very small, which indicates that the inflow from the rivers is less important for the overall water balance of the NCP aquifer system.
Intensive pumping activities have always been blamed to be the main reason for the groundwater depletion in the NCP.They are also held responsible for the problems in the NCP's sustainable water management process.Table 8 and Fig. 11 show that the pumping water takes up about 31 % of the precipitation and irrigation water in the NCP.Therefore, the current pumping amount is not sustainable.The actual pumping taking place is very uncertain and it needs to be addressed  to obtain a more accurate model.Amongst the information which could help improve the model, the pumping rates and the location of wells are the two most important factors.
Figure 13 shows (1) the initial head distribution (1 January 2000) for layer 1, 2 and 3, respectively; (2) the simulated head distribution at the end of the simulation period (31 August 2008) for layer 1, 2 and 3, respectively.For the shallow geological layer (layer 1 in this study area, Fig. 13a  and b), the depression cone near the city Shijiazhuang, which is the capital of Hebei province, had already existed in 2000 and extended further in the next simulation period.While for the deep geological layers (layer 2 and 3, Fig. 13c-f), along with the socio-economic development in the NCP, more and more groundwater is pumped from these two layers to meet the needs of the irrigation and the industrial and agricultural production.As discussed above, irrigation is the biggest consumer of groundwater in the NCP, which is the major reason for the groundwater level decline in this area.The groundwater levels declined about 20 to 30 m in the plain area during the simulation period from 2000 to 2008.

Perspectives on sustainable water management
The results of the water balance analysis have significant implications for the sustainable water management in the NCP.Evapotranspiration is the only actually water consumed or depleted of the hydrologic system (Kendy et al., 2007).Based on this fact, it is an effective and practical way to reduce evapotranspiration to release the water shortage pressure in the NCP area, as about 70 % of the water in the NCP leaves the system by way of ET.Water-saving technologies based on this idea can be adopted to reduce evaporation.And we can change the crop rotation type to help reduce the ET.For example, the crop rotation type with less evaporation is a good choice to replace the existing one.
Groundwater pumping for agricultural irrigation, industrial and domestic use is also listed in Table 8.It shows that irrigation groundwater use accounts for about 79 % of the total pumping over the simulation period (from 83 % in 2000 to 77 % in 2008).At the same time, the percentage of urban use groundwater relative to the total pumping water has increased from approximately 17 % in 2000 to about 25 % in 2007 (23 % in 2008), which means that the industry has developed very quickly and required more water supplies during the simulation period.The industrial development and urbanization in the NCP aggravates the pressure of water supply, as most cities' water supply mainly relies on groundwater.It is obvious that the more emphasis is placed on the economic development, the more water is required and the more pressure will be added to the sustainable development of water resources.
Considering that the water supply in the NCP is mainly groundwater, the protection of groundwater in this area appears especially important.Measures, such as importing surface water from outside river basins and using sources of water other than local groundwater, are effective groundwater protection possibilities.These measures will be greatly beneficial for mitigating the water shortage pressure in the NCP area.Another possibility is to choose measures that can reduce the evapotranspiration, such as changing the crop rotation type or adopting water-saving technologies that reduce evaporation.As shown in Table 8, agriculture is the largest water consumer in the NCP.Therefore it is a primary focus of water conservation attempts to improve the agricultural water use efficiency.Water-saving irrigation techniques, such as border irrigation, low pressure pipe irrigation and sprinkler irrigation (Blanke et al., 2007), are the most effective measures that are used to replace the flood irrigation widely used in the NCP area.
Additional water allocated from the South-to-North Water Diversion Project (SNWDP) will play an important role in solving the water shortage problems in the NCP (Qin et al., 2012).The groundwater pumping for urban use and agricultural irrigation will be reduced once the SNWDP begins to work, as more surface water will be supplied for use.Other than water quantity protection, the protection of the water quality is also very important for the reason that water pollution will aggravate water shortage problems (Qin et al., 2012).The treated sewage can be reused in urban areas and for industry as water supply.Therefore, the wastewater reuse rate should be raised to help alleviate the water deficit of the NCP.A combination of reducing expenditure and controlled wastewater pollution has a dramatic impact on the sustainable water management.This aspect has been examined by Qin et al. (2012) and will be part of the future work of the NCP model.

Perspectives for future work
A key achievement of the present study has been to establish a more reliable water balance for the entire NCP area as a basis for water policy analyses.As other important aspects such as climate change impacts, consideration of water allocation costs, and water management decisions have not been taken into account, the present study should be seen as a platform for further modeling analyses.Scenario analysis is acknowledged as a useful way to utilize hydrological models to support water policy analyses considering alternative assumptions on e.g.future water use, water allocation and climate.Therefore, further work that focuses on climate change, water allocation cost and scenario analysis has already been initiated.

Conclusions
An integrated and distributed hydrological model based on the coupled MIKE SHE/MIKE 11 code has been developed for the entire North China Plain.All major hydrological processes in the land phase of the hydrological cycle such as overland flow, rivers and lakes, unsaturated flow, evapotranspiration, and saturated flow, as well as groundwater pumping and irrigation, are included in the integrated hydrological model.
The model developed in this study is comprehensive and contains several innovative features.The pumping for domestic (urban and rural) and industrial water use and the pumping for irrigation are estimated separately, where the domestic and industrial water consumption is approximated using a system dynamics model incorporating multiple sources of data (Qin et al., 2012).Since groundwater pumping accounts for a substantial fraction in the overall water balance, this improvement is considered significant.Another novel aspect of this study is the comparison of actual ET between model simulation and RS-based estimates as a way of model validation.Moreover the hydrological model used in the present study is an integrated model including both surface water and groundwater components, while similar studies carried out in the same area only consider the groundwater component (Cao et al., 2013) or simply treat the surface water component without considering rivers in the model (Shu et al., 2012).The model is calibrated and validated against time series of hydraulic head data.It is known that groundwater models that are evaluated against groundwater head data alone often may not be able to provide reliable estimates beyond the scope of model calibration.However the potential strength of an integrated model as compared to a groundwater model such as the models developed by Jia et al. (2002), Mao et al. (2005), Wang (2006), Xu (2006), Liu et al. (2008), Wang et al. (2008), Shao et al. (2009), Zhang et al. (2009) and Cao et al. (2013) is that it is possible to use additional data sources other than hydraulic head data to evaluate and/or constrain the model.We compare simulated ET with remote sensing-based ET data (Wu et al., 2012).The model is able to simulate the actual ET in the right order of magnitude and with realistic temporal dynamics, suggesting that the water balance elements simulated by the model are reasonably accurate.Although it has not been possible to validate the stream flow simulations due to the lack of observational data, it is not supposed to have significant impact due to the small quantity compared to the other water balance elements.The results from the comparison of the model-simulated ET and the remote sensing-derived ET are very encouraging and there appears to be good potential in using such data to further constrain the model calibration.For example, Immerzeel and Droogers (2008) presented an innovative approach which incorporates remote sensingderived ET in the calibration of the Soil and Water Assessment Tool (SWAT) in a catchment of the Krishna Basin in southern India.
According to the model-calculated water balance, precipitation, evapotranspiration, irrigation and pumping are the major components controlling the hydrological system inflow and outflow of the NCP.Most of the water (70 %) entering the model area leaves the hydrological system through evapotranspiration, and 29 % leaves the system through pumping.The evapotranspiration is slightly larger than the precipitation, while the groundwater pumping takes up about 31 % of the precipitation and irrigation water in the NCP.Altogether, more outflows than inflows caused the aquifer storage to deplete continuously, which reveals unsustainable water resource development.About 79 % of the pumping water is used as the irrigation water, while the urban use groundwater percentage of the total pumping has increased due to the economic development and population growth.
This model provides a powerful tool for evaluating the water resources in the NCP and particularly how changes in land use and agricultural management impact the groundwater resources.Sustainable water management is a key issue in the process of economic and social development of the NCP.The water balance analysis has significant implications for the sustainable water management in the NCP.Suggestions related to this issue include reducing evapotranspiration, the South-to-North Water Diversion Project, water-saving irrigation techniques and water quality protection.These suggested approaches should be the focus of future studies based on the integrated model developed in this study.Evapotranspiration is the most important outflow component in the overall water balance, and is the only water actually consumed from the hydrologic system.Therefore, effective methods, such as changing crops to less consumptive water users and rainfed crops, fallowing land, and urbanization might be adopted to reduce ET.This will help greatly in mitigating the water deficit pressure and ensuring a sustainable water future for the NCP.In the meantime, measures that used to address water shortage situation in the NCP should consider climate changes and water allocation cost in the future.
) Locations of meteorological stations and rvation wells: the wells marked with rectangles are the selected wells used to show the ration and validation results of groundwater heads; (2) Boundaries of geomorphologic s: (I) piedmont plain, (II) alluvial fan, (III) alluvial plain, and (Ⅳ) coastal plain (Wu et 996; Cao et al., 2013); (3) Locations of capital cities (Beijing, Tianjin, Shijiazhuang Jinan) and some other cities that are included and considered in the NCP model.

Fig. 1 .
Fig. 1.Location and topography of the North China Plain, spanning from 32 • N to 40 • N and from 114 • E to 121 • E(Cao et al., 2013).(1)Locations of meteorological stations and observation wells: the wells marked with rectangles are the selected wells used to show the calibration and validation results of groundwater heads; (2) boundaries of geomorphologic zones: (I) piedmont plain, (II) alluvial fan, (III) alluvial plain, and (IV) coastal plain(Wu et al., 1996;Cao et al., 2013); (3) locations of capital cities (Beijing, Tianjin, Shijiazhuang and Jinan) and some other cities that are included and considered in the NCP model.

Fig. 3 .
Fig.3.Hydrogeological cross section of the NCP along A-A′in Fig.1(adopted from Chen et 996

Fig. 3 .
Fig. 3. Hydrogeological cross section of the NCP along A-A in Fig. 1 (adopted fromChen et al., 2005, andCao et al., 2013).Aquifer units I and II make up layer 1 and also shallow aquifers in the model, while aquifer units III and IV make up layer 2 and 3, respectively; aquifer III and IV are deep aquifers here.

Fig. 4 .
Fig.4.River network of the NCP (IWHR), including 5 major rivers: Luan River, Yongding River, Hutuo River, Zhang River and Wei River.The x and y axis represent the horizontal and vertical directions in unit of meter.

Fig. 4 .
Fig. 4. River network of the NCP (IWHR), including 5 major rivers: Luan River, Yongding River, Hutuo River, Zhang River and Wei River.The x and y axis represent the horizontal and vertical directions in unit of meter.

Fig. 6 .
Fig.6.Scatter plot of averaged simulated and observed groundwater heads for the entire 226 1021

Fig. 6 .
Fig. 6.Scatter plot of averaged simulated and observed groundwater heads for the entire 226 observation wells.Each point represents a value of average simulated/observed head for a well.The points that lack observed head data have been rejected from the figure.The performance criteria values (R, ME, MAE and RMSE) are also shown in the figure.

Fig. 7 .
Fig.7.Monthly comparison of aggregated mean ET values for the entire model domain.RMSE and MAE values given are based on linear regression analysis.

Fig. 7 .
Fig. 7. Monthly comparison of aggregated mean ET values for the entire model domain.RMSE and MAE values given are based on linear regression analysis.

Fig. 8 .Fig. 8 .
Fig.8.Annually aggregated ET from the model and from ET-Watch.The spatial resolution of 1036

Fig. 9 .
Fig.9.Annual mean ET from the model simulation (bold blue line) and the ET-watch (bold red line).The lighter lines represent +/-one standard deviation from the mean values for each of the eight land cover classes considered.

Fig. 9 .
Fig. 9. Annual mean ET from the model simulation (bold blue line) and the ETWatch (bold red line).The lighter lines represent ± one standard deviation from the mean values for each of the eight land cover classes considered.

Fig. 10 .Fig. 10 .
Fig.10.Simulated discharge at one representative location along the five rivers included in 1051

Fig. 11 .
Fig.11.Graphical water balance result of the entire NCP aquifer system.The numbers in the figure are normalized annual averaged values in the type of storage depth with unit of mm/year.

Fig. 11 .
Fig. 11.Graphical water balance result of the entire NCP aquifer system.The numbers in the figure are normalized annual averaged values in the type of storage depth with unit of mm yr −1 .

Table 1 .
(Shu, 2010)the data used in the coupled MIKE SHE/MIKE 11 hydrological model of the NCP(Shu, 2010).

Table 2 .
Shu, 2010)nd elevation information of the meteorological stations (Chinese Meteorological Data Service;Shu, 2010).The average annual reference ET and temperature of each station are also listed here.

Table 4 .
(Shu et al., 2012)ferent soil types(Shu et al., 2012), θ sat is water content at saturation, θ fc is water content at field capacity, θ wp is water content as wilting point and K s is saturated hydraulic conductivity.
hydrogeological unit the associated hydrogeological properties are specified separately.The aquifer is discretized using 27 hydrogeological units (Fig.
is adopted in our study and the hydrogeological parameters are distributed in the layers of the entire NCP.The hydrogeological model is specified via hydrogeological units which are distributed to cover the different layers.For each Hydrol.Earth Syst.Sci., 17, 3759-3778, 2013 www.hydrol-earth-syst-sci.net/17/3759/2013/

Table 6 .
Parameters selected for auto-calibration as well as their final values after the calibration process.K i is horizontal conductivity for geological unit i (i is a positive integer), S y (i) is storage coefficient for geological unit i.

Table 7 .
Performance results of the calibration and validation period.The summary result of mean values comparison shown in Fig.6is also listed here.RMSE: root mean square error; ME: mean error; MAE: mean absolute error; STD res : standard deviation of residuals; and R: coefficient of correlation.