An integrated soil-crop system model for water and nitrogen management in North China

An integrated model WHCNS (soil Water Heat Carbon Nitrogen Simulator) was developed to assess water and nitrogen (N) management in North China. It included five main modules: soil water, soil temperature, soil carbon (C), soil N, and crop growth. The model integrated some features of several widely used crop and soil models, and some modifications were made in order to apply the WHCNS model under the complex conditions of intensive cropping systems in North China. The WHCNS model was evaluated using an open access dataset from the European International Conference on Modeling Soil Water and N Dynamics. WHCNS gave better estimations of soil water and N dynamics, dry matter accumulation and N uptake than 14 other models. The model was tested against data from four experimental sites in North China under various soil, crop, climate, and management practices. Simulated soil water content, soil nitrate concentrations, crop dry matter, leaf area index and grain yields all agreed well with measured values. This study indicates that the WHCNS model can be used to analyze and evaluate the effects of various field management practices on crop yield, fate of N, and water and N use efficiencies in North China.

Conventional field experiments play an important role in assessing the effects of single or multiple factors on crop yield and the fate of N. The processes of soil water dynamics and N cycling in soil-crop systems are complex due to variation in soil properties, weather, and environmental conditions. Models have been successfully used to analyze these complex systems under spatial and temporal variability. Crop and soil models have been used to optimize water and N management practices by simulating water and N dynamics, organic matter turnover and crop growth. Examples of models used around the world include WOFOST 13 , DAISY 14 , HERMES 15 , EPIC 16 , DNDC 17 , RZWQM 18 , DSSAT 19 , APSIM 20 , WNMM 21 , SPACSYS 22 , HYDRUS1D 23 .  and yield (f). K s , saturated hydraulic conductivity; θ s , saturated water content; θ r , residual water content; α, the inverse of the air-entry value; n, pore size distribution index; V n * , maximum nitrification rate; K n , half saturation constant; K d , an empirical proportionality factor; α d Scientific RepoRts | 6:25755 | DOI: 10.1038/srep25755 an equal crop absorption ratio of NH + 4 -N and NO − 3 -N. Each N transformation process was computed as a sink-source term in the CDE, and each of the processes are described detail in next two sections. The boundary conditions (Cauchy type) for the solute (NH + 4 -N and NO − 3 -N) transport equation was used to describe the solute flux at the upper or lower boundary. This CDE was solved by the general upwind difference method, and this procedure effectively avoids numerical dispersion and oscillation even under the conditions of dramatic changes in solute concentration without using dense nodes 39 . Surface broadcast and deep fertilizer applications are regarded as uniformly incorporated within the top 1 cm of soil or at the prescribed application depth (usually 5-10 cm) in the soil, respectively.   Soil organic carbon and N transformation. The module to simulate organic matter turnover was taken from the DAISY model 14,40 , which was originally used for quantitatively evaluating the effect of animal manure on soil-crop systems. The organic matter in soil is divided into three main pools, i.e. dead native soil organic matter (SOM), soil microbial biomass (SMB), and added organic matter (AOM). Each of these distinct pools were considered to contain a continuum of substrate qualities, but to facilitate the description of all turnover processes by first-order kinetic, each of these main pools were divided into two subpools: one with a slow turnover rate (e.g. SOM1, SMB1, and AOM1), and the other with a fast turnover rate (e.g. SOM2, SMB2, and AOM2). The decay rate was proportional to the size of the pool: where ζ i is decomposition or decay rate of pool i (kg C cm −3 d −1 ), C i is carbon content in soil of pool i (kg C cm −3 ), and k i is decomposition or decay rate of pool i (d −1 ). The decomposition or decay rate at actual condition (k AOM ) was derived as the rate at standard abiotic conditions multiplied by abiotic functions taking into account the effects of soil temperature, soil water content, and clay content of the soil. For pools of added organic matter (AOM1 and AOM2), the decomposition rate was calculated from Eq. (2): AOM AOM m m where ⁎ k AOM is the decomposition rate coefficient of AOM at standard conditions (d −1 ), F m (T) and F m (h) are temperature and pressure potential functions, respectively, as detailed by Hansen et al. 14 . For pools of dead native SOM (SOM1 and SOM2), the decomposition rate was computed by

SOM SOM m m m
where ⁎ k SOM is decomposition rate coefficient of SOM at standard conditions (d −1 ), F m (Clay) is clay content function 41 . For pools of microbial biomass (SMB1 and SMB2), the decay rate at the standard conditions as specified in relation to Eq. (3) was assumed to include a specific death rate constant and a specific maintenance rate coefficient: , and ⁎ m is the maintenance coefficient for SMB at standard condition (d −1 ). AOM can be transformed into SMB, and then SMB can be transformed into SOM, which can also be utilized by microorganisms and transformed into SMB. Hansen et al. 40 assumed that soil greenhouse gas emissions and N mineralization-immobilization is closely related to soil microbial activity. Production of CO 2 results from all C-fluxes into the microbial biomass (SMB) pools (substrate utilization efficiencies being less than unity). We assumed a linear relationship between potential denitrification rate and CO 2 emission 42 . Net mineralization rate of N was calculated by the C/N ratio: where S min is net mineralization of SOM (ug cm −3 d −1 ), (C/N) i is C/N ratio of pool i (AOM1, SOM and SMB).
Inorganic nitrogen transformation. Urea hydrolysis. The urea hydrolysis process was computed by a first-order reaction kinetic equation 21 : where N urea is the urea-N concentration in soil (ug cm −3 ), WFPS is the fraction of water-filled pore space (− ), and K urea is the first order kinetic rate constant (d −1 ). Typically urea hydrolysis is completed in several days under hot conditions, while under cold conditions, it takes longer. The models NLEAP, GLEAM, EPIC and others assume urea hydrolysis occurs instantly, and even treated the urea as ammonium directly.
Ammonia volatilization. Ammonia volatilization was simulated using a method proposed by Freney et al. 43  where pH is the value of pH, N am is the concentration of ammonium N in soil (ug cm −3 ), F v (T) and F v (z) are soil temperature functions (°C) and soil depth functions (cm) proposed by Freney et al. 43 .
Nitrification. Nitrification was simulated by the Michaelis-Menten kinetic equation, and modified by soil moisture and temperature 14,40 : Scientific RepoRts | 6:25755 | DOI: 10.1038/srep25755 where V n * is the maximum nitrification rate at 10 °C under optimal soil water condition (ug cm −3 d −1 ), K n is a half saturation constant (ug cm −3 ), F n (T) is the soil temperature function and F n (h) is the pressure potential function proposed by Flowers and O'Callaghan 44 and Addiscott 45 , respectively.
Denitrification. According to Lind 42 , the potential denitrification rate (the extreme anoxic and ample nitrate supplement condition) of the soil can be expressed as a linear function of the CO 2 evolution rate. The actual denitrification rate is determined either by the transport of nitrate to the anaerobic micro sites or the actual microbial activity at these sites. The increased tortuosity when the soil dries leads to discontinuity and thus, denitrification is very limited in dry soil. In the case of ample supply of nitrate, the actual denitrification rate was determined by multiplying the potential denitrification rate by a modified function. Hence, the actual denitrification process was simulated as: where Sden is rate of denitrification (ug cm-3 d-1), α d * is an empirical constant (default value 0.1 g Gas-N/g CO 2 -C); S CO2 is derived from the organic matter model as the evolution of CO 2 from the decomposition of organic matter (ug cm −3 d −1 ); K d is an empirical proportionality factor; N ni is the concentration of nitrate N in soil (ug cm −3 ); F d (θ) is soil water content function described in DAISY model 40 .
Crop growth model. The simulation of crop growth and development stage, LAI, biomass accumulation and allocation, maintenance respiration and growth respiration was computed based on modifications of the PS123 model 46 , which is a generic dynamic crop model to simulate annual crop growth of many crops. The modifications are outlined below.
Crop emergence. The crop relative development stage was divided into two stages: sowing to emergence (stage1) and emergence to maturity (stage2). The first stage was described by a linear function of sowing depth, which was proposed by Mao 47 : = + × Tsum a b depth 1 (11) where Tsum1 is the heat requirement for stage1 (°C), a and b are empirical coefficients, and depth is sowing depth (cm).
Root growth. Root growth and development was computed by Šimůnek et al. 23  where t is time (d), rgr is root growth rate (cm d −1 ), tR min is the initial root growth time (d), tR med , xR med is the time and root depth at the midpoint of development, respectively (cm), xR min and xR max is the initial and maximum rooting depth (cm), respectively, and xR is the root depth (cm) at time t. For overwintering crops, when the temperature in the winter is below zero for 5 continuous days, the roots stop growing. If the soil temperature exceeds 0 °C for 5 continuous days in the spring, the root system begins to grow again.
Root water and nitrogen uptake. Root growth and development was computed by a method proposed by Šimůnek et al. 23 , and root water uptake was calculated by: where T a is actual root water uptake (or crop transpiration) (cm d −1 ), L R is root depth (cm), a w (h, z) is a water stress response function 48 , a s (h ϕ , z) is a salinity stress response function 49 , and b(z) is a root distribution function 48 .
Šimůnek and Hopmans 50 introduced a critical value of water stress index ω c , called the root adaptability factor. This is a threshold value above which root water uptake is reduced in water limited layers of the root zone but the plant compensates by uptaking more water from other layers that have sufficient available water. However, some reduction in potential transpiration occurs below this threshold value, though smaller than that for water uptake without compensation. The water stress index, cf(w), was calculated from Eq. (15).
The shoot and root of crops have different N contents, and the actual root N uptake is determined by the minimum value between the N demand of the crop and soil supply. The crop N stress calibration factor, cf(N), was computed based on the CERES model 38 by: where T s is the accumulated effective temperature from crop emergence (°C); C ANC , C NNC and C MNC are the actual crop N content (%), critical crop N content (%) and minimum crop N content (%), respectively.
The characteristics of WHCNS model. The main characteristics of the integrated model includes: • The primary modules were adapted from existing widely used soil-crop system models, including soil water movement and soil heat transfer routines from HYDRUS1D model 23 , and C and N cycle routines from the DAISY model 14,40 . The crop growth process were based on the PS123 generic crop model 46 , and the gross photosynthetic product was modified by water and N stress calibration factors. The N stress effect on crop growth was adopted from the CERES model. • The numerical convergence problem of the general Richard's equation occurs when heavy rainfall or intensive irrigation (very common in North China) happen, so a simple Green-Ampt model was used to simulate soil water infiltration and water redistribution was simulated using the Richard's equation in the WHCNS model. In addition, to simulate root water uptake, we added a compensatory absorption mechanism to shift water uptake from dry soil zones to wetter soil zones 50 . • The Richard's equation is solved by the Crank-Nicolson implicit method. The convection-dispersion equations of solutes are solved by the general upwind difference method, and this technique can effectively avoid the numerical divergence and oscillation even under the conditions of dramatic changes of solute concentration and without dense nodes 39 . The initial minimum time step is set at 0.1 d, the maximum is set at 1 d, and the number of maximum iteration is set five times. • The WHCNS model can simulate complex and intensive agricultural production systems characteristic to North China, which is typically characterized by conservation tillage, double cropping system, high planting density, film mulching, and intensive water and fertilizer inputs and other management factors. The model allows the user to study the eight irrigation methods defined by FAO56 37 (precipitation, sprinkler, basin, border, every furrow irrigation with narrow/wide bed, alternated furrows irrigation and trickle irrigation) and the model provides four options for fertilizer application (surface, deep, mix and fertigation).

Plot Soil Layer (cm) BD(g cm −3 )
Particle fraction (%)  Table 3. Soil physical and hydraulic properties for soil profiles of the three experiments in Muncheberg 53 . Note: BD is bulk density; θ r is the residual water content; θ s is the saturated water content; α is the inverse of the air-entry value; n is a pore size distribution index; K s is the saturated hydraulic conductivity (l = 0.5).
Scientific RepoRts | 6:25755 | DOI: 10.1038/srep25755 • All input and output files are stored in a user-friendly spreadsheet format, so the simulation results can be directly compared with measured data. The parameters for soil water retention curves can be input as the fitted parameters of the van Genuchten model 49 or as water holding capacity (field capacity and wilting point) 51 . The initial concentration of soil mineral N can be input using a format of mass (mg kg −1 ) or volume (mg L −1 ). • The model was programmed in the C+ + object-oriented programming language, which will allow new features to be added in the future, such as scenario analysis and a parameter optimization module. Recently, a PEST (Parameter ESTimation) module was added to the model to facilitate the optimization of parameters 52 . In addition, a user friendly interface has been designed with the C# programming language within the Microsoft Visual Studio 2008 SDK.
Model Calibration and Validation. The model was tested using the common data sets for comparing models presented in detail by Mirschel et al. 53 collected in Muncheberg, Germany and datasets collected in the North China Plain. The data from Germany were obtained from a 6-year (1992-1998) field experiment carried out at the Centre for Agricultural Landscape and Land Use Research (ZALF) experiment station at Müncheberg, located about 40 km east of Berlin, Germany. There were three rainfed treatments with different fertilizer management in each treatment. Treatment 1 (T1) was intensively managed using inorganic fertilizer and chemical plant protection at a high level; Treatment 2 (T2) was organically managed using only organic manure and non-chemical plant protection and treatment 3 (T3) was an extensive management using a mixture of organic and inorganic fertilizers and chemical plant protection at a low level. The field management and plant and soil N measurement methods were given in detail by Mirschel et al. 53 . The model was calibrated using the treatment 3 datasets, and validated using the treatment 1 and 2 datasets. Calibration consisted of adjusting parameters including crop growth and N transformation in the WHCNS model by comparing the simulated and measured values of soil water content, soil nitrate concentration, crop dry matter (DM) and N uptake during the period from 3 rd Sept. 1992 to 27 th July 1998 in the T3 treatment.
The performance of the WHCNS model was compared to 14 different models (Table 1) which were tested on the Müncheberg data set at the European Conference 54 . Each of these models simulate soil water dynamic, soil N and in some cases, soil C cycles and crop growth using different theories. Some models (Expert-N, SWAP) are designed as toolkits and the users have the choice between different simulation approaches for soil water movement (balance or dynamic method), evapotranspiration and crop growth. For example, three options of crop models, CERES, SPASS, and SUCROS, were linked to Expert-N model, forming three new models denoted as ExN-CER, ExN-SPA and ExN-SUC, respectively. Two models (FASSET and CANDY) combine the capacity approach with soil pore space fractionation for simulating the movement of mobile and immobile water. Most of the models include crop modules, but model approaches range from empirical functions (NDICEA) and simplified temperature-driven approaches (SWIM, WASMOD, CANDY) to more complex mechanical models including photosynthesis, biomass partitioning and root growth development. Crop growth is represented in a generic way in some models (SWAP, Expert-N, HERMES, STAMINA, AGROTOOL, FASSET), but the others use submodels to describe specific crops (CERES, AGROSIM). Depending on their capabilities to simulate crop growth in multiple years, these models were run continuously through the whole crop rotation or were started separately for every application year (CERES, AGROSIM). Eleven models included nitrogen cycle modules, for example, AGROSIM uses a simple N balance approach and zero-order mineralization kinetics, while HERMES describes net mineralization of N from two pools using first-order kinetics. To simulate net immobilization, some models use simple C/N ratios and added organic substances (CERES, SWIM), while the others simulate C and N turnover under more complex conditions. Soil moisture and temperature are the main driving factors for C and N cycles. The denitrification process is included in all of the models except for AGROSIM, STAMINA and AGROTOOL. More detail information about the 14 soil-crop models can be found in book of 'Modelling water and nutrient dynamics in soil-crop systems' 54 .  In order to test the suitability of the WHCNS model in NCP, data collected from four sites (Alxa, Alxa Left Banner in Mongolia; DBW, Dongbeiwang in Beijing, Dawenkou in Shandong province (DWK), and Quzhou in Hebei province (QZ) representing different crops (winter wheat and summer maize), soil types, climate conditions, cropping systems, and field management practices in the NCP (Table 2) were collected or compiled from the literature. The detail experimental design, field management practices and data source in each site are all listed in Table 2.
The Alxa experimental site is located in Left Banner, western Inner Mongolia, China (37.4°-41.9°N, 103.4°-106.9°E), with an elevation of 1150 m. The soils are alluvial mixed with grey desert soils. The region is classified as a warm-temperate desert arid zone with a continental climate. The average annual precipitation is 116 mm, of which 70-80% occurs during the summer season (June to September). Total annual potential evaporation reaches approximately 3005 mm, which is 20 times greater than the annual precipitation. This cropping system is a single crop of maize or wheat produced annually from the middle of April to early October. These crops comprise 80% of the cropped area. The irrigation amounts typically range from 800-1700 mm, and is typically applied through flood irrigation. Typical N fertilizer application rates are approximately 250 kg N ha −1 . The datasets consists of three years (2005,2008,2009) with different water and N management for spring maize 55 .
The DBW experimental station is located in Beijing (39.5°N, 116.2°E), with an elevation of 50 m. The climate is warm-temperate continental monsoon, with an average annual air temperature of 11.5 °C. The average annual rainfall is about 560 mm with 70-80% falling from June to September. The soil is classified as a Cambisol.Surface irrigation by flooding between check banks is widely practiced in the region. The traditional cropping system is winter wheat-summer maize grown within one year. A two-year field experiment under various water and N management was conducted at this site 28 .
The DWK experimental site was located at the Shandong Agricultural University experimental station (36.0°N, 117.1°E) in Tai'an, Shandong Province, China. It is a high productivity region where the average annual grain yield for double cropping systems can approach 15,000 kg ha −1 . The mean annual temperature is 13 °C, with the highest temperature of 26.4 °C in July and the lowest temperature of − 2.6 °C in January. The mean annual rainfall was 697 mm, with the majority rainfall occurring from July to September. A winter wheat and summer maize rotation was used in the experiment, and the soil type was alluvial Cambisols. The field experiment was conducted between October 2009 and September 2012 with different tillage, cultivation, water and N treatments 34 .
The QZ experimental site was located at the China Agricultural University experimental station in Quzhou county, Hebei Province, China (36.6°N, 114.8°E). The county has a continental monsoonal climate. The annual mean air temperature is 13.1 °C. The average precipitation is 556 mm per year, and 70% of the precipitation occurs  between July and September. Total potential evaporation is 1835 mm per year, three times more than annual precipitation. The soil type is alluvial Cambisols and irrigation is provided primarily from groundwater. A two-year field experiment with different water and N treatments was conducted from 1998-2000 56 .
Data from four treatments (Alxa_05_T1, DBW-04-05-T1, DWK-09-10-T1 and QZ-98-99 ) at four sites (Alxa, DBW, DWK and QZ ), respectively were chosen to calibrate the WHCNS model ( Table 2). The soil hydraulic parameters for the model were taken from the literature ( Table 2). Soil C and N transformation parameters came from earlier research at those sites 28,33,[55][56] . The crop parameters were then adjusted according to the measured crop dry weight, LAI, and crop yield by the 'trial-and-error' method. The data from other treatments and years at the four sites were used to validate the model ( Table 2).
Meteorological data, including daily maximum air temperature, minimum air temperature, air humidity, solar radiation, wind speed, and precipitation, were all collected from local weather stations.

Estimating Model Inputs for Calibration.
In order to run the model, inputs need to be estimated. We estimated the inputs for both the European dataset and the NCP dataset using the procedures described below. Model inputs were either estimated or calibrated for the calibration datasets, and then used to test the model performance using the validation datasets.
Soil hydraulic and solute transport parameters. Soil water retention, θ(h), and unsaturated hydraulic conductivity, k(h) functions were estimated by the van Genuchten 49 and Mualem 57 methods, respectively. The parameters for these functions are shown in Table 3. The hydrodynamic dispersion coefficient (D sh , cm 2 d −1 ) in the liquid phase is given by 58 : sh L 0 where D 0 is the molecular diffusion coefficient in free water (cm 2 d −1 ), | | q is the absolute value of the Darcian fluid flux density (cm d −1 ), and D L is the longitudinal dispersivity (cm). Based on the literature 54 , the value of D L was Soil carbon and nitrogen transformation parameters. The microbial activity in soil (soil depth) can be set in the WHCNS model, and the soil N transformation processes (mineralization-immobilization, nitrification and denitrification) are simulated based on it. In this study, the soil microbial activity was set for the top 0-30 cm in the soil. The initial parameters of N transformation were adopted from the default values of the DAISY model, and were all subsequently calibrated according to measured data. The validation parameters were obtained as following:   The long-term dynamic simulation of SOM was not the main goal of this study, and the parameters for SOM decomposition and decay rate were all taken from the earlier studies 14,60 . The initial C/N ratio of residue and distribution coefficient of SOM pools are all from the literature 40 .
Crop parameters. The basic crop parameters were taken from Driessen and Konijon 46 . The partition coefficients and maximum photosynthetic (AMAX) were both calibrated by comparing simulated with the measured total dry weight and LAI. The calibrated crop parameters for the European dataset are shown in Table 4.

Results and Discussion
Sensitivity analysis. The common datasets from Germany were used to analyze the sensitivity of all input parameters of WHCNS. The test parameters can be classified into three categories: soil hydraulic parameters (Table 3), crop parameters (Table 4), and N transformation parameters (listed in section "Soil carbon and nitrogen transformation parameters. "). A sensitivity analysis was carried out by running WHCNS with the value of a single parameter altered by ±10%, while all other selected parameters and variables were constant. The affected output variables included average soil water content, soil nitrate concentration (NO 3 − -N), ammonia concentration (NH 4 + -N) in 0.9 m-soil profile, LAI, DM and crop yield. The results indicated that soil water content was more sensitive to soil hydraulic parameters than crop parameters, and the N transformation parameters had a little effect on soil water content (Fig. 2a). Among soil hydraulic parameters, θ s and n significantly affected soil water content. Sun et al. 61 found that the N transformation parameters typically had less impact on soil water content compared to soil hydraulic parameters and crop parameters. Li et al. 21 also found θ s and n had a high influence on soil water content. Our results are in agreement with these findings.

Soil NO 3
− -N content was more sensitive to crop parameters and soil hydraulic parameters(only n), and the N transformation parameters had little or no effect on NO 3 − -N (Fig. 2b). Soil NH 4 + -N content was more sensitive to N transformation parameters and soil hydraulic parameters, and the crop parameters had a little inflence on NH 4 + -N (Fig. 2c). A change in the n value of van Genuchten equation, V n * and K n significantly affected soil NH 4 + -N content by over 5% (Fig. 2c), indicating that the content of NH 4 + -N was significantly affected by soil pore size distribution index and soil nitrification process. Crop parameters strongly affected LAI and DM, and soil hydraulic parameters and the N transformation parameters had little or no inflence on them (Fig. 2d,e). However, soil hydraulic parameters also strongly affected crop yield, but N transformation parameters had little effect on crop yield (Fig. 2f). Among the crop parameters, T sum , SLA max and AMAX had the highest impact on yield. The main reason was that the value of T sum determines the crop develompent stage, while SLA max is directly related to the LAI. AMAX controls the accumulation of dry matter, and also had a closer relationship with crop yield formation. Richter et al. 62 analyzed the sensitivity of crop parameters for the "Wageningen School" crop model (STAMINA), and found that the most sensitive parameters for yield were related to crop development (T sum in the study) and light interception (SLA max and AMAX). Performance of WHCNS model at the Müncheberg site. We only show the comparison results of simulated and measured soil water contents and nitrate concentrations for because of space limitations. The simulated soil water content at different soil layers in general agreed well with the measured values in T1 (Fig. 3). The comparisons of simulated and measured soil nitrate concentrations in T1 are shown in Fig. 4. The high consistency indicated the ability of the WHCNS model to simulate N transport. Figure 5 shows that the simulated peaks of ammonia (0-30 cm) agreed well with the measured values. The simulated dry matter in three treatments are shown in Fig. 6. The simulation results of all treatments were satisfactory with the exception of the yield of sugar beet in 1993, which was lower compared to measured data. The low simulated dry matter of sugar beet led to a low crop N uptake in 1993 (Fig. 7), and other simulated values matched well with the observed data. The validation results in T1 and T2 indicated that the WHCNS model could be used to simulate the water movement, fate of N and crop growth for different crop rotation systems in Germany.
The performance indices (correlation coefficient, ME, RMSE, IA and NSE ) for the WHCNS model for soil water content and nitrate concentrations at different soil layers in three plots are summarized in Table 5. There was a significant correlation between the simulated and measured. The correlation coefficient ranged from 0.331-0.694. For the soil water contents in the rooting zone (0-90 cm), the ME values ranged from − 0.018-0.001 cm 3 Table 6. The RMSE values of soil water content, mineral N, dry matter and crop N uptake of the WHCNS model were relatively low among 14 compared models, and IA and NSE was higher than most of the other models. Soil mineral N and crop N uptake simulated by WHCNS model had the relatively high NSE value in all simulations (0.38 and 0.79, respectively) compared with the other models. The IA values of dry matter, crop N uptake and soil water content of 0.96, 0.93 and 0.87, respectively, also indicated good model performance. Overall, the WHCNS model gave good agreement with measured data for simulating the dynamics of soil water content and nitrate concentrations in soil profile, as well as crop dry matter and crop N uptake.
There are many reasons that lead to the difference in simulations of soil water movement and N transport. The different approaches for simulating soil water balance (balance or dynamic method) led to differences in simulated soil water content. Krobel et al. 63 compared the performances of the DNDC (balance method) and DAISY (Richard's equation) models in simulating soil water movement and found that the DAISY model had higher accuracy than the DNDC model in simulating soil water content. Crop water uptake is also influenced by root dynamics. Models that simulate more detailed root water uptake generally perform better in simulating soil water content than those with lesser detail 50,64 . Because the compensated root water method was adopted in the WHCHS model, it performed better than most of the other models in simulating soil water content. For N transport, the complexity of N transformations considered in the model affects the simulated result of soil mineral N content. Eleven of the models included N cycle process modules 54 . AGROSIM uses a simple N balance approach and zero-order mineralization kinetics, while HERMES describes net mineralization of N from two pools using first-order kinetics. For net immobilization, some models use simple C/N ratios and added organic substances (CERES, SWIM), while others simulate C and N turnover under more complex conditions (such as soil moisture and temperature as the main driving factors for C and N cycles) 54 . In additional, those influencing factors of soil water movement and crop growth will also affect N transformation and transport. In this study, the integrated model (WHCNS) combined the best aspects of widely used soil-crop models, and improved the ability of model to represent various environmental conditions, and gave better performance compared to most of the other 14 soil-crop system models. Figure 11 shows the validation results of soil water content, nitrate concentrations, dry matter and LAI in four sites in the NCP. The coefficients of determination between simulated and measured soil water content and nitrate concentrations are in the range of 0.42-0.87 and 0.42-0.71, respectively at four sites. For the validation of soil water content, the RMSE values ranged from 0.03-0.04 cm 3 cm −3 (Table 7); the IA indices were all over 0.78; and the NSE values had the lowest value of 0.34 in Alxa, and values of 0.78, 0.76 and 0.63 in DWK, QZ and DBW, respectively. Hu et al. 29 applied the RZWQM2 model to simulate soil water content for double cropping system in NCP, and the RMSE range was 0.03-0.07 cm 3 cm −3 . Van Liew and Garbrecht 65 analyzed several soil water models and suggested that the simulation performed well when the value of NSE > 0.36 and IA > 0.7. In this study, these statistical indices are all within these reported acceptable ranges.  The validation results for the nitrate simulation in different soil layers showed that all the IA indexes were more than 0.77, the RMSE values ranged from 2.58-7.04 mg kg −1 , the ME indices ranged between − 1.23 and − 0.28 mg kg −1 , which showed the model underestimated the soil nitrate content. The NSE values ranged from 0.29-0.63 for the four sites. This might be due to the complex N transformation process. Kersebaum et al. 54 compared the simulation results of 13 soil-crop models, and found that the dynamic of mineral N usually had a negative NSE value (− 0.81-− 0.20). Hu et al. 29 reported that the value of NSE in general was negative, which was caused by the high variability observed for soil NO 3 − -N measurement. Given the complexity of N transformation, the ranges of these indices appear to be acceptable. These results indicated that the model performed well in simulating soil water dynamic and nitrate transport under different water and N management practices.

Calibration and evaluation of WHCNS in North China. Soil water and nitrogen.
Crop growth. The species of crop growth simulated at the four sites included winter wheat, summer maize and spring maize. The validation results showed that the model explained more than 84% of the variation in LAI and dry matter (Fig. 11), and all the slopes of the regression lines were close to 1, except for the LAI in Alxa (1.5) due to limited measured data (only four measurements), while the RMSE values for dry matter ranged from 1006-2962 kg ha −1 (Table 7). Total dry matter had the largest error in QZ (2962 kg ha −1 ), and the smallest error in DWK (1006 kg ha −1 ). Figure 12 shows the simulation results of crop yields at the four sites for the validation datasets. In general, the simulation results matched well with observed values. The determination coefficients for the four sites was 0.94 (Fig. 12). The RMSE values ranged from 243-1097 kg ha −1 , and the IA indices were all over 0.70 (Table 7). Palosuo et al. 66 compared the performance of eight widely used soil-crop system models for winter wheat yield simulation at seven sites in Northern and Central Europe, and found that the DAISY and DSSAT models gave better performance, the values of RMSE were the lowest (1428 and 1603 kg ha −1 , respectively) and IA (0.71 and 0.74, respectively) was the highest. Among those models, CROPSYST underestimated crop yields (ME, − 1186 kg ha −1 ), while HERMES, STICS and WOFOFST overestimated crop yields with ME values of 1174, 1272 and 1213 kg ha −1 , respectively. In another comparison of nine crop models for spring barley yield simulation at seven sites in Northern and Central Europe, Rotter et al. 67 found that HERMES, MONICA and WOFOST outperformed other models in simulating crop yield with the lowest RMSE values of 1124, 1282 and 1325 (kg ha −1 ), respectively. Compared to those indices simulated by these models, the WHCNS performed reasonably well in simulating crop yield in North China. Figure 13 shows the box-plot for the simulated and measured crop yields of winter wheat, summer maize and spring maize. The results showed that the data distribution of the simulated total and specific crop yields was similar to that of the measured ones but with less variance. This issue was also found by many researchers when simulating the crop growth at a regional scale [68][69] . This was probably due to lower sensitivity of the model to the N fertilizer application because there are large amount of accumulated mineral N in soil profile in North China 28,34,56,70 . Another reason may be the uncertainty of the observed yield caused by the extreme climate, diseases and pests [66][67] . In our study, the extreme climate led to an abnormally low yield of winter wheat in 2011 (3220 kg ha −1 ) in DBW. However, the current models do not include modules to simulate impacts of these factors and further study is needed.
The WHCNS model successfully simulated soil water content, nitrate concentrations, crop yields, and LAI in the North China, and explained the difference in crop yield under different water and N management practices. Hence, the model has the potential to be used in simulating soil water movement, N cycle, and crop development under different climate conditions, soils, crop rotations and a variety of water and fertilizer management practices in NCP.
However, as a newly developed model, more testing is needed using long-term datasets from different sites to assess the uncertainty caused by climate change. In terms of functions, the new model cannot simulate the effects of other nutrients (such as phosphorus fertilizer), diseases and pests on crop growth. In addition, the greenhouse gas (CH 4 ), dissolved carbon (DOC) and the dissolved N (DON) play an important role in lowland agricultural production. The WHCNS model cannot currently simulate these processes. In order to simulate hydrologic processes on a regional scale, a hydrological and groundwater module should be incorporated in the model in the future.

Conclusion
The open access dataset from a field experiment in Germany, available after a modeling workshop in Europe was used to test the newly developed WHCNS model. We then compared the results with the simulation results of 14 other soil-crop system models. The WHCNS model ranked in the top three based on the NSE and IA indices for soil water contents (NSE, 0.36; IA, 0.87), soil nitrate concentrations (NSE, 0.38; IA, 0.79), crop dry matter (NSE, 0.84; IA, 0.96) and crop N uptake (NSE, 0.36; IA, 0.87). We concluded that WHCNS gave good performance in simulating soil water dynamic, nitrate and ammonia transport, crop growth development and grain yield under different crop rotation and fertilizer management practices.
In addition, the WHCNS model was also tested using datasets from four different sites in North China across different soils, climate conditions, cropping systems, and intensive water and fertilizer management practices. The IA and NSE values for LAI and dry matter were close to 1. The determination coefficient for the crop yields ranged from 0.84-0.99. And the RMSE values ranged from 243-1097 kg ha −1 . The IA indices were all over 0.70. All these results indicated that the newly developed WHCNS can be used to analyze and evaluate the grain yield, fate of N, WUE and NUE in the intensively cultivated farmland in North China.