Lihong HEXingguo MO*Shi HU and Suxia LIU

1 Key Laborator y of Water Cycle and Related Land Sur face Processes，Institute of Geographic Sciences and Natural Resources Research，Chinese Academy of Sciences，Beijing 100101(China)

2 College of Resources and Environment，University of Chinese Academy of Sciences，Beijing 100049(China)

3 Sino-Danish Center，University of Chinese Academy of Sciences，Beijing 100049(China)

ABSTRACT Accurately simulating the soil nitrogen(N)cycle is crucial for assessing food security and resource utilization efficiency.The accuracy of model predictions relies heavily on model parameterization.The sensitivity and uncertainty of the simulations of soil N cycle of winter wheat-summer maize rotation system in the North China Plain(NCP)to the parameters were analyzed.First，the N module in the Vegetation Interface Processes(VIP)model was expanded to capture the dynamics of soil N cycle calibrated with f ield measurements in three ecological stations from 2000 to 2015.Second，the Morris and Sobol′algorithms were adopted to identify the sensitive parameters that impact soil nitrate stock，denitrif ication rate，and ammonia volatilization rate.Finally，the shuffled complex evolution developed at the University of Arizona(SCE-UA)algorithm was used to optimize the selected sensitive parameters to improve prediction accuracy.The results showed that the sensitive parameters related to soil nitrate stock included the potential nitrif ication rate，Michaelis constant，microbial C/N ratio，and slow humus C/N ratio，the sensitive parameters related to denitrif ication rate were the potential denitrif ication rate，Michaelis constant，and N2 O production rate，and the sensitive parameters related to ammonia volatilization rate included the coefficient of ammonia volatilization exchange and potential nitrif ication rate.Based on the optimized parameters，prediction efficiency was notably increased with the highest coefficient of determination being approximately 0.8.Moreover，the average relative interval length at the 95%conf idence level for soil nitrate stock，denitrif ication rate，and ammonia volatilization rate were 11.92，0.008，and 4.26，respectively，and the percentages of coverage of the measured values in the 95%conf idence interval were 68%，86%，and 92%，respectively.By identifying sensitive parameters related to soil N，the expanded VIP model optimized by the SCE-UA algorithm can effectively simulate the dynamics of soil nitrate stock，denitrif ication rate，and ammonia volatilization rate in the NCP.

Key Words: ammonia volatilization，denitrif ication rate，global sensitivity analyses，shuffled complex evolution developed at the University of Arizona(SCE-UA)algorithm，vegetation interface processes model

INTRODUCTION

Increase in grain yield relies heavily on the large investment of nitrogen(N)fertilizer in the North China Plain(NCP)，which is a key region of grain production(Hartmannet al.，2015;Wanget al.，2017).The average annual application of N fertilizer is as high as 550—600 kg ha-1in the winter wheat-summer maize rotation system，which is much higher than the demand of crops in this region(Chen，2003).Additionally，improper management practices，such as excessive irrigation，incorrect application rates，and inadequate sowing rates，will limit the improvement in N use efficiency(NUE)and lead to a series of environmental problems(Zhu and Chen，2002;Juet al.，2006;Liuet al.，2006;Guoet al.，2010).Studying the mechanisms of N transformation in the winter wheat-summer maize rotation system is essential to achieve high crop yield，efficiently utilize nitrogenous resources，and reduce the risk of environmental pollution(Zhouet al.，2006;Erismanet al.，2007;Padgett and Minnich，2008;Zhang Xet al.，2017;Zhanget al.，2018).These are of great signif icance for developing rational agricultural policies to ensure food security and environmental sustainability.

Due to the constraints of complex soils，crops，and planting patterns，f ield experiments are time-consuming and laborious(Smithetal.，2008).Conversely，process-based models，which are important research tools，can economically assess the effects of interactions between soil，climate，and agricultural management on soil N dynamics.To date，numerous models have been developed，including the EPIC(Williams，1995)，APSIM(McCownetal.，1996)，DAISY(Abrahamsen and Hansen，2000)，RZWQM(RZWQM Development Teamet al.，1998)，WNMM(Li Yet al.，2007)，DNDC(Liet al.，1992)，and WOFOST(van Diepenet al.，1989).These models were typically designed for specif ic objectives and applied to multiple types of studies in the NCP(Lu and Fan，2013;Zhang Jet al.，2017;Chen Z Wet al.，2018).For example，the WOFOST is a mechanistic model that forecasts crop production，biomass，or water use efficiency based on biochemical processes(e.g.，photosynthesis and respiration)inf luenced by environmental conditions and is widely used to quantify potential crop yield in the NCP(Wuet al.，2006;Lu and Fan，2013).The DNDC，designed to simulate carbon(C)and N biogeochemistry in agroecosystems，is used to predict soil C dynamics(Zhang Xet al.，2017;Chen Zet al.，2018)，N leaching(Liet al.，2014)，and trace gas emissions(Li Het al.，2010).The RZWQM coupled with the DSSAT crop models is initially focused on assessing the productivity and water quality of various cropping systems for soil，climate，and f ield management(RZWQM Development Teamet al.，1998)，and its applications include N transport and NUE in the NCP(Yuet al.，2006;Liet al.，2015).Due to the focuses and limitations of each model，previously developed models have advantages in assessing crop yield，soil C dynamics，water movement，NUE，and emissions of trace gases to a certain extent.To accurately describe the processes of soil N cycle，it is necessary to develop a soil N cycle module that is suitable for the climatic conditions and f ield management in the NCP.

Generally，mechanistic models require many input parameters，which make it difficult to calibrate and lead to an obvious uncertainty.It is necessary to identify parameters that have a signif icant impact on model performance through sensitivity analysis(SA)before calibration.Fixing non-inf luential parameters at reasonable values will greatly decrease the computational cost(Saltelli and Annoni，2010).Most previous studies have used local sensitivity analysis(LSA)methods to evaluate model responses by consecutively varying one parameter while keeping other parameters at a constant value(Heinen，2006;Moet al.，2012;Lianget al.，2016).For example，Lianget al.(2016)analyzed soil nitrate and ammonium concentrations by varying±10%of each input parameter in the WHCNS model，and found that soil ammonium content was sensitive to N conversion，while the hydraulic parameters had little effect.However，the local SA method cannot analyze the effect of interactions between the different parameters on model performance.However，global sensitivity analysis(GSA)methods，such as the Morris(Morris，1991)，Sobol′(Sobol′，2001)，and FAST(Cukieret al.，1973)algorithms，can be used to evaluate the response to a model performance by simultaneously varying all parameters，which helps to fully understand the sensitivity of model parameters(Lamboniet al.，2009;Moet al.，2012;Lianget al.，2017).The Morris and Sobol′algorithms were used to evaluate the responses of nitrate leaching to soil hydraulic and soil organic matter transformation parameters.It was found that saturated soil water content and f ield capacity were the most inf luential parameters(Lianget al.，2017).Chen Z Wet al.(2018)used the Sobol′method to screen inf luential parameters for N transportation and identif ied potential denitrif ication rate as the most sensitive parameter.However，few studies have adopted the GSA methods to completely identify the sensitive parameters of soil N cycle.

Generally，the sources contributing to modelling uncertainty are incorrect input data，unreasonable model structure，and uncertain parameters values(Refsgaard and Storm，1990).Various methods have been developed to deal with parameter uncertainty，such as the generalized likelihood uncertainty estimation(GLUE)(Beven and Binley，1992)，simulated annealing algorithm(Kirkpatricket al.，1983)，genetic algorithm(Whitley，1994)，and the shuffled complex evolution developed at the University of Arizona(SCE-UA)algorithm(Duanet al.，1994).These algorithms have been widely used for parameter calibration and uncertainty analysis(Vrugtet al.，2003;Daiet al.，2009;Zhanget al.，2009).For example，six parameters of CN-SIM model were carried out with a multi-objective genetic algorithm(i.e.，NSGA-II)，by minimizing the relative root mean squared error(RMSE)between observations and simulations with the data of respired C and soil inorganic N measured for three soil types(Cavalli and Bechini，2012).Raf iqueet al.(2015)used PEST software to optimize the inf luential parameters for N2O emissions in the DayCent model and found that the simulated accuracy of the optimized N2O f lux was improved by 63%.Previous studies have individually optimized parameters related to the N cycles，such as soil organic matter decomposition or denitrif ication(Cavalli and Bechini，2012;Raf iqueet al.，2015).Based on previous SA results，selecting the appropriate calibration algorithm can greatly improve the efficiency of parameter calibration.However，few studies attempted to optimize parameters related to soil N cycle as a whole in an agroecosystem.In this study，based on SA results，the SCE-UA algorithm was used to calibrate model and provide technical support for the promotion and application of the VIP(Vegetation Interface Processes)model.

The VIP model mainly accounts for soil organic matter decomposition，photosynthesis，soil hydrothermal processes，vegetation dynamics，and energy partition(Mo and Liu，2001;Moet al.，2012).After years of development and uncertainty studies of model parameters，the VIP model has been used successfully to simulate winter wheat-summer maize yield，evapotranspiration，and water use efficiency in the NCP(Moet al.，2012).This study attempted to identify sensitive parameters that affect soil N cycle in a winter wheat-summer maize rotation system in the NCP and evaluate the uncertainty interval of the simulation results caused by parameter uncertainty.The objectives of this study were:i)to expand the N cycle module in the VIP model based on the existing soil N cycle theory and methodology;ii)to conduct a global sensitivity analysis using the Morris and Sobol′algorithms to identify sensitive parameters of the N cycle in the study area;and iii)to optimize parameters using the SCE-UA algorithm based on previous SA results and evaluate the uncertainty interval of the VIP performance.

MATERIAL AND METHODS

Description of the study area

The NCP(31°14′—40°25′N，112°48′—122°45′E)is located in eastern China(Fig.1).The area of the NCP is about 330 000 km2，70%of which is agricultural land.It belongs to the semi-arid monsoon climate zone with most rainfall occurring in summer(Hu，2012).The annual average temperature varies from 10 to 15°C，and the sunshine duration is about 2 800 h in the central region.It is a key granary in China，with a dominant cropping system of a winter wheat-summer maize rotation(Hu，2012).

Fig.1 Map of the study area and spatial distributions of agricultural ecology experimental stations.

Data

Datafor simulation.Daily meteorological records，soil data，and f ield management data were employed(Table I).Among these，considerable data on daily meteorological variables collected at three ecological stations(Fengqiu，Yucheng，and Luancheng)were from the Chinese Ecosystem Research Network(CERN)(http://www.cnern.org.cn/data/).The study period was 2007—2010 at Fengqiu，2000—2015 at Yucheng，and 2003—2010 at Luancheng.Soil data includes soil physical，chemical，and hydraulic properties in different layers(Table II).Soil physical and chemical properties(including particle fraction，bulk density，and organic matter)were collected from the literature(Ma，2004;Li X Xet al.，2007;Huanget al.，2015).The pedo-transfer function was used to estimate soil hydraulic characteristics from the bulk density and the percentages of sand，silt，and clay(Schaapet al.，2001).Indirect estimation of soil hydraulic parameters using the pedo-transfer function has been widely used in the studies of agricultural hydrological simulation(Mishraet al.，2013;Jianget al.，2015;Li and Ren，2019).The planting system was a winter wheat-summer maize rotation.The planting methods and f ield management at the three ecological stations were collected from the literature(Fanget al.，2006;Huang，2011)and CERN.Details on planting date，harvesting date，Nfertilization，and irrigation are presented in Table III.Urea was used as N fertilizer.Nitrogen was divided into base fertilizer and topdressing application.Irrigation was performed immediately after N topdressing both in the wheat and maize seasons.Irrigation was provided if there was no precipitation or the rainfall was low.

TABLE I Overview of site information collected to drive the expanded Vegetation Interface Processes(VIP)model for the simulations of soil N cycle in a winter wheat-summer maize rotation system

Data for calibration and validation.Observed data were measured at three ecological stations(Table IV).Field measurements of soil nitrate stock(kg ha-1)were obtained from CERN.To determine soil nitrate content，10 g of each soil sample was extracted with 50 mL of 2 mol L-1KCl.The nitrate concentration of the extract was measured at 220 and 275 nm(Normanet al.，1985).The nitrous oxide(N2O)f lux was determined in a closed static system.During the measurement，a gas sample was taken from a sampling port at the top of the box，and 20 mL was taken each time.At the same time，the temperature in the box，soil temperature at 10 cm soil depth，and air temperature were measured.The samples were then analyzed using an HP-5890 gas chromatography/electron capture detector(GC/ECD)(Hewlett-Packard，USA).The measurement conditions were as follows:the front column and analysis column were stainless steel packed columns，and the inner diameters were 2 and 1 mm;column temperature was 55°C and ECD temperature was 330°C;and high purity N(30 mL min-1)was used as the carrier gas.The ammonia volatilization f lux was measured using a semi-open static system，as described by Liuet al.(2003).After fertilization，the ammonia capture system was placed in each plot.A vacuum pump was used to make the ammonia volatilized from the soil entering the scrubber with the air stream，and then absorbed by 2%boric acid.Absorbed ammonia was determined by indophenol blue colourimetry at 625 nm(UV1800，Shimadzu，Japan)(Denmeadet al.，1976).

TABLE II Soil physical，chemical，and hydraulic properties a)in different layers at three ecological stations(Fengqiu，Yucheng，and Luancheng)used for the simulations of soil N cycle in a winter wheat-summer maize rotation system by the expanded Vegetation Interface Processes(VIP)model

TABLE III Overview of site information relating to planting methods and f eild management，collected for the simulations of soil N cycle in a winter wheat-summer maize rotation system，at three ecological stations(Fengqiu，Yucheng，and Luancheng)

TABLE IV Observation information collected from three ecological stations(Fengqiu，Yucheng，and Luancheng)to calibrate and validate the expanded vegetation interface processes(VIP)model for the simulations of soil N cycle in a winter wheat-summer maize rotation system

Model description

The processes in the original VIP model(Mo and Liu，2001;Moet al.，2012)include photosynthesis，vegetation growth dynamics，soil hydrothermal movement process，soil C cycle，and energy and water balance.Soil heat transfer and soil water movement are calculated using the thermal diffusion equation and the Richards equation，respectively.Soil Ncycle module accounts for turnover and mineralization between the soil organic Npools with the conversion between inorganic N pools neglected.On this basis，the processes of soil N cycle module including nitrif ication，volatilization，denitrif ication，and immobilization were expanded in this study(Details shown in Fig.S1，see Supplementary Material for Fig.S1).Below are the main processes of C and N cycles of the model.

Soil organic matter mineralization.Soil organic matter mineralization module was adopted from the CENTURY model(Partonetal.，1993).Soil organic C and N are divided into f ive pools:structural litter pool，metabolic litter pool，microbial biomass pool，slow humus pool，and passive humus pool.Different organic C pools are controlled by the f irstorder kinetic equation(Eq.1)and the decomposition rate is proportional to the size of the pools(Partonet al.，1993):

whereSC，iis the C stock of the pooli(kg ha-1);tis the time required to decompose the pool(d);andkiis the decomposition rate of the pooli(d-1)，which is inf luenced by soil moisture and temperature，as shown in Eq.2:

whereis the decomposition rate of organic matter under optimal moisture and temperature conditions(d-1);Tis the soil temperature(°C);wis the soil water content(g kg-1);andf(T)andf(w)represent temperature and moisture correction functions，respectively.The decomposition rates of different organic N pools are shown in Eq.3:

whereSN，iis the soil N stock of pooli(kg ha-1)and(C/N)iis the C/N ratio of pooli.

Nitrif ication rate.The substrate(NH+4，O2，and CO2)content，soil temperature，and water content are considered to be the abiotic factors that affect the nitrif ication rate.In the case of relatively high soil moisture(1.5＜pF＜2.5，pF=log(-h)，wherehis the soil water potential)，pH value，and temperature，soil microbial metabolic activity is limited by soil organic C content.In the absence of oxygen stress，most ammonium-N is rapidly oxidized into nitrate，and the nitrif ication rate is described by the Michaelis-Menten equation(Eq.4)(Abrahamsen and Hansen，2000):

whereδnitis the nitrif ciation rate(μg cm-3d-1);is the nitrif ication rate constant under optimal temperature and moisture content(μg cm-3d-1);Camis the ammonium concentration(μg cm-3);Knitis the semi-saturation constant(μg cm-3);and(T)and(h)are the soil temperature(T)and water potential(h)functions，respectively.

Denitrif ication rate.Denitrif ication rate is a dimensionless multiplicative function of soil potential denitrif ication rate that explicates nitrate concentration，soil temperature，and moisture effects(Johnssonetal.，1991).Hénault and Germon(2000)proposed a NEMIS denitrif ication rate prediction model based on the measured values of denitrif ication rate of undisturbed soil，nitrate concentration，and water-f illed pore space.Denitrif ication rate is expressed as Eq.7(Hénault and Germon，2000):

whereξdenis the denitrif ication rate(μg cm-3d-1);is the denitrif ication rate constant under the optimal temperature and moisture content(μg cm-3d-1);Kdenis the semisaturation constant(μg cm-3);(T)is the soil temperature function;(wf)is the water-f illed pore space function;andCntis the soil nitrate content(μg cm-3).

Ammonia volatilization rate.The rate of ammonia volatilization is affected by soil ammonium concentration，soil temperature，and soil pH(Sherlock and Goh，1984).It is expressed as Eq.8:

whereφvolis the ammonia volatilization rate(kg ha-1d-1);NHxis the total amount of soil NH4and NH3(kg ha-1);Tis the soil temperature(°C);pH is the soil pH value;ksoil-airis the exchange coefficient of ammonia in soil and air;andkaq-soilis the exchange coefficient of ammonia in soil liquid water and soil.

Urea hydrolysis rate.Most of the N fertilizer applied in the NCP is urea.Under high temperature and soil moisture conditions，urea hydrolysis occurs over a few days，while it takes a longer time under lower temperature and dry conditions(Li Yet al.，2007).Models，such as NLEAP，GLEAMS，and EPIC，assume that urea hydrolysis occurs immediately and that urea is treated directly as ammonium-N(Lianget al.，2016).The urea hydrolysis rate(Rhys)is described by Eq.11(Li Yet al.，2007):

whereCureais the urea content in soil(kg ha-1);WFPS is the water-f illed pore space;andKureais the f irst-order kinetic rate constant(d-1).

Model initialization and parameters

Daily meteorological data measured at the three ecological stations，including maximum temperature，minimum temperature，average temperature，relative humidity，sunshine hours，precipitation，and wind speed，were used to drive the VIP model.The initial values of vegetation characteristics，such as initial crop height，maximum rooting depth，soil temperature，and soil water content，were specif ied based on f ield observations or general site knowledge(Mo and Liu，2001;Ma，2004;Li X Xet al.，2007;Moet al.，2012;Huanget al.，2015;CERN).Soil organic matter turnover and mineralization parameters in the expanded VIP model include organic matter decomposition rate，C/N ratio，ammonia volatilization rate，and denitrif ication rate(Table V)，which were obtained from previous studies(Sherlock and Goh，1984;Williamset al.，1989;Partonet al.，1993;Abrahamsen and Hansen，2000;Hénault and Germon，2000;Li Yet al.，2007).The model was spun-up for 300 years，using randomly selected forcing from 5 years of observation data at Yucheng.The C and N pools of the model stabilized after about 200 years.Following spin-up，the model was continually run forward from October 1，2007 to September 30，2010 at Fengqiu，from October 1，2000 to September 30，2015 at Yucheng，and from October 1，2003 to September 30，2010 at Luancheng，with the observed forcing data.Outputs obtained from these simulations were compared with the corresponding observations.

Sensitivity analysis and optimization

For complex agroecosystem models，parameter optimization requires a high computational cost.The parameters，which have considerable effects on model performance，can be identif ied through SA.Then，non-inf luential parameters are f ixed at reasonable values，which will decrease the computational cost without reducing model accuracy(Campolongoet al.，2007).In this study，the Morris method(Morris，1991)was f irst used to screen the inf luential parameters related to soil nitrate stock，denitrif ication rate，and ammonia volatilization rate in the expanded VIP model.Then，the f irst-order，second-order，and total sensitivity indices of the identif ied parameters were analyzed quantitatively using Sobol′method(Sobol′，2001).Finally，the SCE-UA algorithm(Duanet al.，1994)automatically optimized the selected parameters within the range of values，in order to achieve optimal agreement between simulated and corresponding measured values of soil nitrate stock，denitrif ication rate，and ammonia volatilization rate.The software Problem Solving Environment for Uncertainty Analysis and Design Exploration(PSUASE)was used for SA and optimization(Tong，2017).This provides several global SA methods，including the Morris，Sobol′，and SCE-UA algorithms，which were utilized in this study.

TABLE V Soil and crop N parameters in the expanded VIP model for the simulations of soil N cycle in a winter wheat-summer maize rotation system

Objectivefunction.During the optimization process，the SCE-UA algorithm is used to solve the minimum value of the objective function，which is used to evaluate the consistency between the observed and simulated values(Duanet al.，1994).Different objective functions are used to evaluate the different characteristics of the model.Guptaet al.(1998)summarized several objective functions for parameter optimization.In this study，the RMSE was used to def ine the objective function.The observations were grouped and weighted accordingly depending on type.The objective function was divided into three groups:soil nitrate stock，ammonia volatilization rate，and denitrif ication rate.The specif ic form of the objective function was as follows:

whereφis the objective function;nni，nam，andndenare the numbers of observed values of soil nitrate stock，ammonia volatilization rate，and denitrif ication rate，respectively;OandSrepresent the measured and simulated values，respectively;αni，αam，andαdenare the weighting coefficients of soil nitrate stock，ammonia volatilization，and denitrif ication rate，respectively，and their values were 1/3 in this study.

Evaluating indicators of model performance.It is difficult to comprehensively evaluate model performance using single parameters and it is necessary to describe the simulation accuracy by different indicators.In this study，three statistical indicators were used to assess the agreement between measured and corresponding simulated values.The indicators are RMSE(Willmott，1982)，Nash-Sutcliffe eff iciency coefficient(NSE)(Nash and Sutcliffe，1970)，and coefficient of determination(R2)(Nagelkerke，1991).In addition to the statistical indicators(RMSE，NSE，andR2)noted above，two other statistical indicators were used:average relative interval length(ARIL)(Eq.13)and percentage of measured values coverage at the 95%conf idence level(P-95CI)(Eq.14)，as proposed by Li Let al.(2010)to evaluate the uncertainty of the model performance.

whereS97.5%andS2.5%are the quantile values at 97.5%and 2.5%of the simulated value set after 10 000 samples by Monte Carlo andNinis the number of measured values within the 95%conf idence interval.A small ARIL value indicates a very narrow uncertainty interval and a large P-95CI value indicates the reliability of this interval(Li Let al.，2010).

RESULTS

Sensitivity analysis

Screeningsensitiveparametersby Morrismethod.The Morris method is adopted to distinguish qualitatively sensitive parameters related to soil nitrate stock，denitrif ication rate，and ammonia volatilization rate.The value of parameterp(Yang，2011)in the present study was 8.Additionally，this study assumed that all inputs were evenly distributed within the parameter range.The results are presented in a plot with modif ied mean values(Fig.2).Large deviation of a parameter from the origin of coordinates indicates its importance.Sensitive parameters related to soil nitrate stock included potential nitrif ication rate，urea decomposition rate，semi-saturation constant，the decomposition rate of structural litter，the decomposition rate of metabolic litter，microbial decomposition rate，microbial C/N ratio，and slow humus C/N ratio.Sensitive parameters related to denitrif ication rate were potential denitrif ication rate，potential nitrif ication rate，the decomposition rate of structural litter，the decomposition rate of metabolic litter，microbial C/N ratio，slow humus C/N ratio，and the Michaelis constant of nitrif ication.Sensitive parameters related to ammonia volatilization rate were ammonia volatilization exchange coefficient，maximum nitrif ication rate，and urea decomposition rate.

Fig.2 Sensitivity analysis results of Morris method for soil nitrate stock(a—c)，denitrif ication rate(d—f)，and ammonia volatilization rate(g—i)at three ecological stations，Fengqiu(a，d，and g)，Yucheng(b，e，and h)，and Luancheng(c，f，and i).The sample size(N)in this study，calculated by N=(n+1)·r，was 440(n and r are the numbers of parameters and replication，respectively).See Table V for the detailed descriptions of the abbreviations.

Quantitativeanalysisof sensitiveparametersusing Sobol′method.Because it is difficult to determine the quantitative sensitivity of specif ic parameters using the Morris method，the Sobol′method is used to estimate the f irst-order，second-order，and total sensitivity indices of expanded VIP parameters.The f irst-order and total sensitivity indices of soil nitrate stock，denitrif ication rate，and ammonia volatilization rate at Fengqiu，Yucheng，and Luancheng are shown in Fig.S2(see Supplementary Material for Fig.S2).At the three stations，the ranking of parameters for soil nitrate stock based on the results of f irst-order sensitivity index was consistent with that based on the total sensitivity indices.The potential nitrif ication rate was the most sensitive parameter(Fig.S2A)，and its f irst-order sensitivity indices were 0.202，0.246，and 0.189 at Fengqiu，Yucheng，and Luancheng，respectively.This result can be explained by the Michaelis-Menten equation(Eq.4)(Abrahamsen and Hansen，2000)that potential nitrif ication rate determines the nitrif ication process directly.The results of total sensitivity indices(Fig.S2B)indicated that soil nitrate stock was highly affected by potential nitrif ication rate，urea decomposition rate，semi-saturation constant，and metabolic litter decomposition rate，which accounted for average 29%，16%，14%，and 11%of soil nitrate stock variability，respectively.Based on the f irst-order index，the denitrif ication rate was most sensitive to potential denitrif ication and nitrif ication rates at the three ecological stations(Fig.S2C).The f irst-order indices were 0.19，0.25，and 0.19 for potential denitrif ication rate and 0.18，0.24，and 0.17 for potential nitrif ication rate at Fengqiu，Yucheng，and Luancheng，respectively.The total sensitivity index(Fig.S2D)revealed that the denitrif ication rate was sensitive to potential denitrif ication rate，potential nitrif ication rate，and urea decomposition rate，which accounted for average 31%，23%，and 30%of denitrif ication rate variability，respectively.Based on the results of f irstorder sensitivity index，ammonia volatilization exchange coefficient was the most sensitive parameter for ammonia volatilization rate at all stations(Fig.S2E)，and its f irst-order indices were 0.17，0.20，and 0.17 for Fengqiu，Yucheng，and Luancheng，respectively.The sensitive parameter obtained from the total sensitivity indices was similar to that from f irst-order sensitivity index，except one difference:exchange coefficient between soil and soil water(vw2s)rather than exchange coefficient between soil and air(vs2a)was the most sensitive parameter for ammonia volatilization rate at Fengqiu and Yucheng(Fig.S2F).In general，the ammonia volatilization rate was sensitive to ammonia volatilization exchange coefficients(vs2aandvw2s)，potential nitrif ication rate，and urea decomposition rate，which accounted for average 23%，28%，17%，and 12%of ammonia volatilization rate variability，respectively.The second-order index represents the inf luence of the interaction between two parameters on the objective function.Fig.3 shows the second-order index between different parameters related to soil nitrate stock，denitrif ication rate，and ammonia volatilization rate at Fengqiu，Yucheng，and Luancheng.For soil nitrate stock，the second-order index between potential nitrif ication rate and urea decomposition rate demonstrated the most interaction，being 0.004 6，0.002 6，and 0.003 5，at Fengqiu，Yucheng，and Luancheng，respectively.For denitrif ication rate，the interaction between potential denitrif ication rate and urea decomposition rate was the most sensitive，and the second-order indices were 0.014 5，0.011 6，and 0.0015 9，respectively.For ammonia volatilization rate，the interaction between two ammonia volatilization exchange coefficients(vs2aandvw2s)was the most sensitive，and the second-order indices were 0.008 9，0.008 6，and 0.008 4，respectively.When the f irst-order index of a parameter on the objective function was signif icant，the second-order interaction was also signif icant.

Fig.3 Sensitivity analysis results of Sobol′second-order indices for soil nitrate stock(a—c)，denitrif ication rate(d—f)，and ammonia volatilization rate(g—i)at three ecological stations，Fengqiu(a，d，and g)，Yucheng(b，e，and h)，and Luancheng(c，f，and i).See Table V for the detailed descriptions of the abbreviations.

Model calibration by SCE-UA optimization algorithm

Based on the period and source of observed data collected at the three ecological stations(Fengqiu，Yucheng and Luancheng)，soil nitrate stock，denitrif ication rate，and ammonia volatilization rate were calibrated to validate performance of the VIP model in soil N cycle.Additionally，to improve calibration efficiency and accuracy，the SCE-UA algorithm was adopted to optimize the selected parameters based on the prior SA results，so that it achieved a good f it between the measured and corresponding simulated values of soil nitrate stock，denitrif ication rate，and ammonia volatilization rate.Evaluation indicators of soil nitrate stock，denitrif ication rate，and ammonia volatilization rate between reference values and the SCE-UA optimization algorithm are shown in Table VI.Specif ic parameter selection and the range of values for each parameter at the three ecological stations are shown in Table VII.

Dynamics of soil nitrate stock in the winter wheatsummer maize rotation system.Accurate simulation of soil nitrate stock helps to further assess the dynamics of soil organic matter turnover and mineralization.In this study，the simulated soil nitrate stocks were evaluated with f ieldobserved values in the soil surface(0—20 cm)and bottom layer(20—120 cm)during the growth periods of wheat and maize.The datasets of soil nitrate stock from 2007 to 2008 at Fengqiu and from 2000 to 2001 at Yucheng were used to calibrate the expanded VIP model.The datasets from 2008 to 2009 at Fengqiu and from 2001 to 2002 at Yucheng were used to validate the model(Fig.S3，see Supplementary Material for Fig.S3).Based on SCE-UA algorithm，the simulated soil nitrate stocks were close to the observed values in soil surface layer(0—20 cm)(Fig.S3A-1，A-2)，but overrated compared to the measured values in the bottom layer(20—120 cm)(Fig.S3A-3，A-4)at Fengqiu.The values of RMSE，NSE，andR2were 9.54，0.59，and 0.61，respectively.The validation also obtained reasonable statistical results(RMSE=11.52，NSE=0.43，andR2=0.49).At Yucheng，the values of soil nitrate stock were not signif icantly overestimated or underestimated in both calibration(Fig.S3B-1，B-3)and validation processes(Fig.S3B-2，B-4).The VIP model show good performance in simulating soil nitrate stock，with an RMSE value of 8.36，NSE of 0.68，andR2of 0.62 in the calibration，and an RMSE value of 10.76，NSE of 0.59，andR2of 0.53 in validation.Following calibration by SCE-UA algorithm，the RMSE values of soil nitrate stock simulation were reduced by 34.2%and 31.6%at Fengqiu and Yucheng，respectively，and theR2values were increased by 41.2%and 39.1%，respectively.As shown in Fig.S3，the applied urea was rapidly decomposed and converted into nitrate-N，resulting in a rapid increase in soil nitrate stock.The soil nitrate stock decreased slowly during the seeding to reviving stage，and then decreased rapidly after winter wheat reviving.This was mainly because ammonia volatilization and denitrif ication intensif ied with increasing temperature，and crop N uptake increased rapidly after the reviving stage.These results indicated that following calibration by SCEUA algorithm，the VIP model could simulate soil nitrate stock well.Given the complexity of N transformation，the ranges of these indices(RMSE，NSE，andR2)appeared to be acceptable.These results indicated that the model performed well at simulating soil nitrate transport in this region.Lianget al.(2016)reported that the NSE value of soil nitrate stock simulation in different soil layers ranged from 0.29 to 0.63.Kersebaumet al.(2007)compared the simulation results of 13 soil-crop models and found that the NSE value of mineral N was usually negative.In the present study，these indices were all within these reported acceptable ranges.

TABLE VI Evaluation indicators a)of soil nitrate stock，denitrif ication rate，and ammonia volatilization rate between reference values and shuffled complex evolution developed at the University of Arizona(SCE-UA)optimization algorithm

TABLE VII Model-related parameters optimized by shuffled complex evolution developed at the University of Arizona(SCE-UA)algorithm at three ecological stations，Fengqiu，Yucheng，and Luancheng

Denitrif icationrate.The datasets of denitrif ication rate during the growth periods of wheat and maize from 2008 to 2009 at Fengqiu and Luancheng were used to calibrate the VIP model，and datasets from 2009 to 2010 were used to validate the model(Fig.S4，see Supplementary Material for Fig.S4).Based on SCE-UA algorithm，the statistical results of denitrif ication rate showed that the VIP model could effectively simulate denitrif ication rate during the calibration(Fig.S4A-1，B-1)and validation(Fig.S4A-2，B-2).The RMSE，NSE，andR2values were 0.04，0.71，and 0.66 in calibration at Fengqiu.The corresponding values for Luancheng were 0.03，0.77，and 0.72，respectively.The validation also showed fair results based on the RMSE(0.06and 0.04)，NSE(0.62 and 0.64)，andR2(0.53 and 0.68)at Fengqiu and Luancheng，respectively.After calibration by SCE-UA algorithm，at Fengqiu and Luancheng，the RMSE values of denitrif ication rate simulation were reduced by 29.8%and 21.5%，and theR2values were increased by 33.2%and 35.6%，respectively.Besides，the 1:1 line of the f itting diagram also showed that the simulation of denitrif ication rate was acceptable.The probability and intensity of denitrif ication were greater during the maize growing season compared with the wheat growing season，as shown in Fig.S4.The low water content is the main factor limiting the denitrif ication process during the wheat growing period.Furthermore，although the intensity of denitrif ication increased，it did not reach the peak in 2008 at Fengqiu(Fig.S4A-1)，because the concentration of nitrate was relatively low.During the maize growing period in 2009，the intensities of denitrif ication were similar at both Fengqiu(Fig.S4A-2)and Luancheng(Fig.S4B-1，B-2).This was because the substrate(soil nitrate stock)was maintained at a relatively high level and was conducive to microbial denitrif ication.In this study，theR2value of denitrif ication rate simulation ranged from 0.53 to 0.71，which is similar to previous reports(Li Yet al.，2007;Guet al.，2016).The DNDC model was used to simulate N2O emissions during the winter wheat-summer maize growing season and the correlation coefficient ranged from 0.79 to 0.90.Guet al.(2016)evaluated nitrous oxide emissions from a winter wheat-summer corn rotation system during a 25-year fertilization trial in northwestern China.In addition，the VIP model generated a higher denitrif ication rate in the maize season than that in the winter wheat season，which is consistent with previous results.

Ammonia volatilization rate.Ammonia volatilization is one of the main ways of Nfertilizer loss in agroecosystems.The datasets of ammonia volatilization rate during the growth periods from 2008 to 2009 at Fengqiu and from 2003 to 2004 at Luancheng were used to calibrate the VIP model;the datasets from 2009 to 2010 at Fengqiu and Luancheng were used to validate the model.Figure S5(see Supplementary Material for Fig.S5)presents the simulated and measured values of ammonia volatilization rate under fertilization treatments based on the SCE-UA algorithm and the reference values at Fengqiu and Luancheng.At Fengqiu，based on SCE-UA algorithm，the simulated ammonia volatilization rate was close to the measured values in model calibration，with an RMSE value of 0.91(Fig.S5A-1);however，it was overestimated during the model validation process，with an RMSE value of 1.77(Fig.S5A-2).The NSE values were 0.59 and 0.43 in model calibration and validation，respectively，and theR2values were 0.64 and 0.55，respectively.At Luancheng，the predicted ammonia volatilization rates were close to those observed in both the calibration(Fig.S5B-1)and validation(Fig.S5B-2)processes.A good simulation of the ammonia volatilization rate was achieved with an RMSE value of 0.86，NSE of 0.68，andR2of 0.79 in model calibration.The validation also showed good results(RMSE=0.98，NSE=0.59，andR2=0.73).The RMSE values of ammonia volatilization rate simulation were reduced by 31.5%and 30.2%，respectively，at Fengqiu and Luancheng and theR2values were increased by 35.3%and 36.2%，respectively，following optimization by SCE-UA algorithm.The results showed that the ammonia volatilization rate was lowest in mid-January and highest in early August during the wheat and maize growing seasons.Since relatively high temperature and humidity are benef icial to urea hydrolyzation，the volatilization rate gradually decreases and tends to stabilize several days after application of N fertilizer.Besides，N fertilizer is generally applied to soil surface，making it easily volatilized.Modeled ammonia f luxes are highly sensitive to fertilization practices.As indicated in Fig.S5，fertilizer application virtually had impact on ammonia f luxes.Cuiet al.(2014)modeled trace gas f luxes after fertilizer application and showed that fertilizer application enhanced ammonia f luxes.The inf luence of fertilizer application date could be linked to two factors:change in temperature with time and the lag between fertilizer application and planting dates.Fertilization before planting increased ammonia f luxes，because the time for nitrif ication and N loss to the atmosphere is increased.

Uncertainty intervals of simulation per formance

To quantitatively evaluate simulation uncertainty caused by parameter uncertainty，10 000 parameter sets were sampled with the Monte Carlo method.The ARIL and P-95CI of soil nitrate stock，denitrif ication rate，and ammonia volatilization rate are used to evaluate the uncertainty of parameters.The predicted uncertainties of soil nitrate stock，denitrif ication rate，and ammonia volatilization rates at Fengqiu are presented in Fig.4.The results showed that the uncertainty ranges of the simulation results included most observations.The ARILs at the 95%conf idence level for soil nitrate stock，denitrif ication rate，and ammonia volatilization rate were 11.92，0.008，and 4.26，respectively.Approximately 68%，86%，and 92%of observed values for soil nitrate stock，denitrif ication rate，and ammonia volatilization rates were within the 95%conf idence interval in the VIP model.This indicated that individual parameter uncertainty may be more important than other sources of uncertainty.Parameter uncertainty accounted for most of the deviation between the observed and the simulated values.

Fig.4 Uncertainty intervals of soil nitrate stock(S nit)at depths of 0—20 and 20—120 cm，denitrif ication rate(ξden)，and ammonia volatilization rate(φvol)at the Fengqiu Ecological Station from Oct.1，2008 to Jun.1，2009(winter wheat growing period)and from Jun.1，2009 to Oct.1，2009(summer maize growing period).There were 10 000 parameter sets sampled by the Monte Carlo method.Arrow represents fertilization practice.SCE-UA=shuffled complex evolution developed at the University of Arizona.

DISCUSSION

Effectiveness of sensitivity analysis

Sensitivity analysis is a key step in the identif ication of sensitive parameters.It is also an effective way to elucidate how the parameters affect model performance and provides a scientif ic basis for further parameter optimization(Vrugtet al.，2003).In this study，the Morris screening method was used to screen the parameters that inf luence simulations of soil nitrate stock，denitrif ication rate，and ammonia volatilization rate.The theory of Morris is based on the fundamental effect，which represents the model response to a particular parameter.The effect of each parameter is evaluated by the mean of the fundamental effects of each parameter(Campolongoet al.，2007).

The sensitivity results obtained by the Morris method indicated that nitrif ication and decomposition of soil organic matter are the main sources of soil nitrate.Similarly，some studies have shown that nitrif ication and soil organic matter decomposition are the main processes affecting soil N concentration(Xenakiset al.，2008;Lianget al.，2016).The denitrif ication rate is sensitive to the potential denitrif ication and nitrif ication rates.The Michaelis-Menten equation(Eq.7)revealed that the potential denitrif ication rate and soil nitrate stock are two important factors affecting the denitrif ication process，and the potential denitrif ication rate is the main parameter affecting the soil nitrate stock，which indirectly conf irms the results of sensitivity analysis.According to the results of the sensitive parameters of ammonia volatilization，it can be known that ammonia volatilization exchange coefficient，maximum nitrif ication rate，and urea hydrolysis rate have great inf luences on ammonia volatilization process.In the ammonia volatilization control equation，the gaseous ammonia concentration in the atmosphere was assumed to be zero;therefore，ammonia volatilization exchange coefficient will signif icantly affect ammonia volatilization rate.These analysis results demonstrate that the Morris method can effectively screen sensitive parameters.

In general，the inabilities to quantitatively estimate the contribution of parameters to the output variance and distinguish the interactions between different parameters are the main limitations of the Morris algorithm(Yang，2011).The second objective of SA is to quantitatively evaluate the inf luence of the selected sensitive parameters on model performance.The Sobol′method was adopted to comprehensively evaluate the sensitive parameters selected by the Morris method.Sobol′sensitivity indices are analysis of variance-like methods，which are adopted to quantitatively assess the total indices of the parameters.In the present study，the sample sizes of the Morris and Sobol′methods were 440 and 2 300，respectively，and the computational cost of the Sobol′method was about f ive-times that of the Morris method.The SA results in this study indicated that the potential nitrif ication rate had the greatest inf luence on soil nitrate stock.Lianget al.(2016)evaluated the SA of the WHCNS model in the NCP and showed that N conversion parameters，such as maximum nitrif ication rate and half-saturation constant，were important factors affecting soil nitrate stock.Notably，the results of the SA are obtained at a specif ic time scale and parameter range and depend on the parameterization scheme of the analysis method(Confalonieriet al.，2010).

Effectiveness of SCE-UA algorithm on parameter optimization

Few studies have investigated the optimization of parameters related to the N cycle of agroecosystems model by SCE-UA algorithm，which is used widely in hydrological models(Gelleszunet al.，2017;Kanet al.，2018).Optimization of selected parameters by SCE-UA algorithm signif icantly increased the simulation accuracy of nitrate stock，denitrif ication rate，and ammonia volatilization rate.The NSE for soil nitrate stock simulation increased to 0.67 after parameter optimization in this study.By comparing the performances of multiple models，it was found that the NSE value of soil nitrate stock simulation varied from-0.81 to 0.20(Kersebaumet al.，2007).The VIP model simulation has a relatively higher NSE value.For the denitrif ication rate simulation，following parameter optimization by SCE-UA algorithm in this study，the NSE andR2values increased to 0.73 and 0.78，respectively，which was similar to previous reports(Li Yet al.，2007).The correlation coefficient of N2O emissions varied from 0.60 to 0.71 during the winter wheat-summer maize growing season using the DNDC model(Li Hetal.，2010).These results show that the VIP model has a higher simulation accuracy than DNDC model.For the ammonia volatilization rate simulation，theR2value was 0.62 after parameter optimization by SCE-UA algorithm in this study.Dubacheetal.(2019)used the DNDC95 model to simulate the dynamics of ammonia volatilization after urea application;theR2value varied from 0.44 to 0.59，which is lower than the results of the present study.This indicates that，after optimized by SCE-UA algorithm，the VIP model can effectively capture the patterns in the seasonal ammonia volatilization during the wheat and maize growing periods.

Model parameter values

In this study，the potential nitrif ication rate and its halfsaturation constant，with the range of 6.4—9.6 and 44—66 g m-3d-1，respectively，had higher sensitivity for soil nitrate stock simulation，and the optimized values were 9.5 and 44.0 g m-3d-1，respectively.Lianget al.(2016)applied the WHCNS model to agroecosystems in the NCP and optimized the N conversion parameters by PEST.They found that the potential nitrif ication rate and its half-saturation constant were 10.2 and 51.7 g m-3d-1，respectively，which were slightly higher than the values obtained in the present study.Parameter optimization is closely related to the range of parameter values.The ranges of the above two parameters in the WHCNS model are 10—30 and 25—75 g m-3d-1，respectively，which are larger than those in this study during the optimization process.This may explain why the optimized values of the above two parameters are larger than those in this study.Furthermore，the objective function and the difference in the optimization algorithm also explain the different optimized values.The denitrif ication rate simulation in the VIP model is based on the NEMIS scheme(Hénault and Germon，2000)，in which the potential denitrif ication rate and its half-saturation constant are important parameters inf luencing the denitrif ication rate，and the optimized values by SCE-UA algorithm were 2.46 and 22.3 kg ha-1d-1，respectively.Guet al.(2016)simulated the denitrif ication rate of a maize-wheat rotation system in southwestern China using the NOEv2model and showed that the denitrif ication rate was 2.96 kg ha-1d-1，which is similar to the value of this study.From the results of parameter optimization，it was found that the optimized values of some parameters were unchanged.Further research is needed to determine how to set the threshold of parameter sensitivity and determine the range of optimization parameter values.

CONCLUSIONS

In this study，Morris and Sobol′algorithms were used to perform a global sensitivity analysis of the N cycle module.The selected sensitive parameters，based on the global sensitivity analysis results，were optimized by SCE-UA algorithm to improve the simulation accuracy of expanded VIP model.The results of sensitivity indices indicated that Morris and Sobol′algorithms could effectively identify sensitive parameters of soil N cycle module at three ecological stations.The sensitive parameters related to soil nitrate stock were potential nitrif ication rate，Michaelis constant，microbial C/N ratio，and slow humus C/N ratio.The sensitive parameters related to denitrif ication rate were potential denitrif ication rate，Michaelis constant，and N2O production rate.The sensitive parameters for ammonia volatilization rate included ammonia volatilization exchange coefficient and potential nitrif ication rate.After optimization by SCE-UA algorithm，the VIP model with the expanded N module can effectively simulate soil nitrate stock，denitrif ication rate，and ammonia volatilization rate during the winter wheat and summer maize growing periods.Furthermore，the optimized value of each parameter was within a reasonable range and more than 68%of observed values were within the 95%conf idence interval.

ACKNOWLEDGEMENTS

This work was f inancially supported by the National Natural Science Foundations of China(Nos.41790424 and 41471026).We are grateful to CERN for providing soil N data.We also thank the editors and reviewers for their assistance with this manuscript.

SUPPLEMENTARY MATERIAL

Supplementary material for this article can be found in the online version.