Using a modified DNDC biogeochemical model to optimize field management of a multi-crop (cotton, wheat, and maize) system: a site-scale case study in northern China

It is still a severe challenge to optimize the field management practices for a multi-crop system when simultaneously aiming at yield sustainability and minimum negative impacts on climate as well as atmosphere and water quality. This site-scale case study was devoted to developing a biogeochemical process model-based approach as a solution to this challenge. The best management practices (BMPs) of a three-crop system growing cotton and winter wheat–summer maize (W–M) in rotation, which is widely adopted in northern China, were identified. The BMPs referred to the management alternatives with the lowest negative impact potentials (NIPs) among the scenarios satisfying all given constraints. The independent variables used to determine the NIPs and those utilized as constrained criteria were simulated by the DeNitrification-DeComposition model, which was modified in this study. Due to the unsatisfactory performance of the model in daily simulations of nitric oxide (NO) emission and net ecosystem exchange of carbon dioxide (NEE), the model was modified to (i) newly parameterize the soil moisture effects on NO production during nitrification, and (ii) replace the original NEE calculation approach with an algorithm based on gross primary production. Validation of the modified model showed statistically meaningful agreements between the simulations and observations in the cotton and W–M fields. Three BMP alternatives with overlapping uncertainties of simulated NIPs were screened from 6000 management scenarios randomly generated by Latin hypercube sampling. All of these BMP alternatives adopted the baseline (currently applied) practices of crop rotation (3 consecutive years of cotton rotating with 3 years of W–M in each 6-year cycle), the fraction of crop residue incorporation (100 %), and deep tillage (30 cm) for cotton. At the same time, these BMP alternatives would use 18 % less fertilizer nitrogen and sprinkle or flood-irrigate ∼23 % less water than the baseline while adopting reduced tillage (5 cm) for W–M. Compared with the baseline practices, these BMP alternatives could simultaneously sustain crop yields, annually enlarge the soil organic carbon stock by 4 ‰ or more, mitigate the aggregate emission of greenhouse gases, NO release, ammonia volatilization, and nitrate leaching by ∼7 %, ∼25 %, ∼2 %, and ∼43 %, respectively, despite a ∼5 % increase in N₂O emission. However, further study is still necessary for field confirmation of these BMP alternatives. Nevertheless, this case study proposed a practical approach to optimize multi-crop system management to simultaneously achieve multiple United Nations Sustainable Development Goals.

1 Introduction

Globally, fiber crops (i.e., cotton) and cereals such as wheat and maize have long played a relevant role in human society because they are primary sources of materials for the textile and food industries. In China, although cotton cultivation only covers 2.0 %–3.9 % of the annual crop harvest area (cotton lint production totalled 5.3 × 10⁶–7.6 × 10⁶ Mg from 2007 to 2016), the cultivation of cereals is quite large. Wheat and maize accounted for 39 % and 26 % of the harvest area and represented 1.29 × 10⁸ and 2.2 × 10⁸ Mg of grain, respectively, in 2016 (China Statistical Yearbook, 2017).

Northern China is both the second most important area of cotton production and the largest region of the winter wheat–summer maize double-cropping system (i.e., both crops are harvested within a year, and they are hereinafter referred to as W–M) in the country (e.g., Cui et al., 2014). Crop rotations of cotton and W–M have commonly been grown in this region, alternating every 3 to 5 years (e.g., Liu et al., 2010, 2014). Over the last few decades, the yields of cotton, wheat, and maize have been increased by employing intensified agricultural management practices, such as increased fertilizer inputs, advanced irrigation methods, and so on (e.g., Han, 2010). A recent study indicated that the cotton cropping system in northern China persistently functioned as an intensive carbon or net aggregate greenhouse gas (GHG) source compared with W–M (Liu et al., 2019). These previous studies have revealed that the change in the storage of soil organic carbon (ΔSOC), net ecosystem aggregate GHG emissions (NEGE), and other biogeochemical processes involving the multiple cropping system in northern China are likely closely related to the rotation pattern of cotton and W–M (e.g., Liu et al., 2010, 2014, 2019; Lv et al., 2014).

To maintain high productivity, the three-crop rotation system for cotton and W–M in northern China are characterized by large additions of synthetic nitrogen fertilizers and irrigation amount (e.g., Chen et al., 2014; Galloway et al., 2004), at 60–140 and 550–600 kg N ha⁻¹ yr⁻¹, and 140–200 and 90–690 mm yr⁻¹ for cotton and W–M, respectively (e.g., Ju et al., 2009; Liu et al., 2014; Wang et al., 2008). High nitrogen and water inputs can result in high release potentials for nitrogenous pollutants, which can induce a series of environmental problems, such as water pollution due to increased nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) leaching (e.g., Collins et al., 2016). In addition, other field management practices, e.g., tillage and crop residue treatment, can also affect the emissions of reactive nitrogen and contribute to negative environmental effects (e.g., Zhang et al., 2017; Zhao et al., 2016). Therefore, the evaluation of the multiple-cropping system (e.g., rotations of cotton and W–M) is shifting from a single-goal method aimed at increasing crop yields to a multi-goal approach (e.g., Cui et al., 2014; Garnett et al., 2013; Zhang et al., 2018). This multi-goal strategy aims to simultaneously sustain/increase crop productivity to ensure food security, increase SOC content to improve soil fertility, mitigate NEGE to alleviate climate warming, reduce ammonia (NH₃) volatilization and nitric oxide (NO) emission to secure air quality, and abate ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching to protect water quality.

According to the multi-goal approach, an objective method is applied to identify the best management practice (BMP), which evaluates each decision variable (using price-based proxies or other measures) and screens the best option with the minimal negative impact potential (NIP) under the given constraints at the annual scale (e.g., Cui et al., 2014; Xu et al., 2017). To screen the BMP, it is essential to quantify the biogeochemical effects of management practices at the annual scale. Field experiments are often only capable of focusing on the decision variables of very few management practices during short periods (e.g., Ding et al., 2007; Liu et al., 2010, 2015; Wang et al., 2013a, b). However, this limitation of field experiments can potentially be overcome by process-oriented biogeochemical models, such as DeNitrification-DeComposition (DNDC) (Li et al., 1992; Li, 2000, 2007, 2016), DAYCENT (Delgrosso et al., 2005), and LandscapeDNDC (Haas et al., 2012).

A three-crop (cotton, winter wheat, and summer maize) system in southern Shanxi Province was selected for this model-based site-scale case study. This study was undertaken to (i) diagnose problems with the DNDC95 model version that was validated in Cui et al. (2014) against the comprehensive field measurements of the selected W–M fields; (ii) make modifications and then validate the modified model for both cotton and W–M cropping systems, especially for the variables that determine NIPs and those involved in given constraints; and (iii) investigate the biogeochemical effects of various rotation patterns with different field management practices and identify the multi-goal BMP alternatives based on the modified model simulations. These efforts were undertaken to test two hypotheses: that a validated process-oriented biogeochemical model is capable of addressing a challenging issue – the optimization field management practices of a three-crop rotation system; and that the field management of an intensive three-crop rotation system can be optimized to simultaneously sustain the current crop yields and annually increase the SOC stock by 4 ‰ so as to implement the international “4 Per 1000” initiative (https://www.4p1000.org/, last access: 31 May 2019) – an action of the Paris Agreement aimed at combating climate change, mitigating aggregate greenhouse gas emissions, and reducing other negative impacts on the environment.

2 Materials and methods

2.1 Brief introduction to DNDC95

The original DNDC95 model used by Cui et al. (2014) and in this study is one of the latest DNDC versions (http://www.dndc.sr.unh.edu/model/GuideDNDC95.pdf, last access: 15 September 2017). This model consists of two components with six modules in total. Driven by the given primary ecological factors, the first component simulates the field states of a soil–plant system, such as the soil chemical and physical status, vegetation growth, and organic matter decomposition. Driven by the soil-regulating variables yielded by the first component, the second component simulates the core biogeochemical processes of carbon and nitrogen transformations and the physical processes of liquid and gas transport; thus, the annual dynamics of net ecosystem exchanges of carbon dioxide (CO₂)(NEE); emissions of methane (CH₄), nitrous oxide (N₂O), NH₃, and NO; and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching and the inter-annual dynamics of SOC and NEGE are simulated. These features enable the model to investigate the integrative biogeochemical effects of the rotation patterns of multiple crops and/or other management practices based on comprehensive validation. The minimum inputs used to facilitate the model simulation include (i) the meteorological variables of daily precipitation and maximum/minimum temperature; (ii) the soil (cultivated horizon) properties of the clay fraction, bulk density, SOC content, and pH; (iii) the crop parameters for the yield potential, thermal degree days (TDD) for maturity, and the mass fractions and carbon-to-nitrogen (C∕N) ratios of the grain, root, and leaf plus stem; (iv) the management practice variables of sowing and harvest (dates), fraction of incorporated/retained residue at harvest, tillage (date and depth), irrigation (date, method, and water amount), and fertilization (date, type, method, nitrogen amount, and C∕N ratio of organic manure); and (v) other variables (annual means of the NH₃ concentration in the atmosphere and the ammonium plus ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration in rain water). For more details regarding the model, please see Li et al. (1992) and Li (2000, 2007, 2016).

2.2 Problem diagnosis and model modification

The daily simulations of the original DNDC95 showed poor model performance for NO emissions from the cotton field, e.g., with a very negative Nash–Sutcliffe efficiency index (NSI) of −1.03 for the 333 observations in 2009. Meanwhile, the original model often failed to capture the daily high NEE fluxes on rainy or cloudy days despite the agreement between the simulation and observation of the annual cumulative NEE (Cui et al., 2014). According to the program codes, a constant fraction (0.003) of a nitrification rate (F_n; in g N m⁻² d⁻¹) was adopted in the original model to calculate the daily NO production in the nitrification process. This was found to account for the former problem as the fraction could vary with soil moisture, which was mechanically similar to the N₂O production in nitrification (shown by the model program codes). The latter problem was ascribed to the algorithm adopted to calculate daily NEE. In the original model, a daily NEE flux was calculated as the residue of daily CO₂ release by soil heterotrophic respiration and daily CO₂ uptake by the increase in cumulative net primary production (NPP). The daily cumulative NPP was calculated by portioning the input of potential crop yields to the day following a given NPP growth curve (shown by the model program codes; Li, 2016). Consequently, the insensitivity of a daily NPP increase to radiation intensity reduction resulted in a more negative daily NEE on a rainy or cloudy day. In order to solve these two problems, the model was modified in this study as follows.

In the model version modified in this study, the effect of the soil moisture (SM) in water-filled pore space (WFPS; dimensionless 0–1 fraction) on NO production was parameterized by referring to that for N₂O production during nitrification and incorporated into the function by replacing the constant fraction mentioned above (Eq. 1). This modification was adopted to reflect that high soil moisture facilitates the production of NO as a by-product of nitrification (NO_p; in g N m⁻² d⁻¹). The maximum fraction of NO production (K_n; dimensionless 0–1 fraction) during nitrification was calibrated as 0.03 using the observed daily NO fluxes from October 2007 to October 2008 in the cotton field. The other observations of daily and annual NO fluxes from both adjacent lands with different field treatments (Table S1 in the Supplement) were used to validate this modification.

$$\begin{array}{}\text{(1)}& {\mathrm{NO}}_{\mathrm{p}}={\mathrm{SM}}^{\mathrm{5.0}}{K}_{\mathrm{n}}{F}_{\mathrm{n}}\end{array}$$

In the model version modified in this study, a daily NEE flux (g C m⁻² d⁻¹) is calculated as the residue of the daily CO₂ release from ecosystem respiration (ER; in g C m⁻² d⁻¹) and daily CO₂ uptake due to gross primary production (GPP; in g C m⁻² d⁻¹). While the daily ER is simulated as it is in the original model, the daily GPP is calculated using Eq. (2), which is a widely applied hyperbola function of photosynthetically active radiation (PAR; in mmol m⁻² d⁻¹; e.g., Wang et al., 2013a; Zheng et al., 2008).

$$\begin{array}{}\text{(2)}& \mathrm{GPP}=\mathit{\alpha }\mathrm{PAR}\cdot {\mathrm{GPP}}_{max}/(\mathit{\alpha }\mathrm{PAR}+{\mathrm{GPP}}_{max})\end{array}$$

The photosynthetic parameters in Eq. (2), α (g C mmol⁻¹) and GPP_max (g C m⁻² d⁻¹), denote apparent quantum yield and maximum GPP, respectively. Each parameter is quantified as the product of shoot standing biomass (B_s; in g C m⁻²) and biomass-specific apparent quantum yield (f₁; in g C mmol⁻¹ (g C m⁻²)⁻¹) or specific GPP_max (f₂; in g C m⁻² d⁻¹ (g C m⁻²)⁻¹) corrected by an adaptation factor (a; dimensionless) reflecting inter-annual variations of a crop (Eqs. 3, 4). The variable B_s is simulated as it is in the original model. The seasonal dynamics of f₁ and f₂ is a function of the normalized plant developing stage (ds; dimensionless 0–1 fraction) (Eqs. 5, 6). The functions of f₁ and f₂ take the forms presented by Zheng et al. (2008) for winter wheat, whereas their empirical parameters of a, b, c, d₁, d₂, g, h, i, j, l, and m can be calibrated to adapt both functions to given conditions. Two daily NEE fluxes per week were randomly selected from the year-round eddy covariance observations in both cotton and W–M cropping systems (Cui et al., 2014; Liu et al., 2019; Wang et al., 2013a) to calibrate the values of these parameters specifically for cotton, winter wheat, and summer maize, while remaining daily data (independent NEE observations) were used to evaluate this modified NEE algorithm. For the calibrated values of these parameters, refer to Table S2.

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathit{\alpha }=a{f}_{\mathrm{1}}{B}_{\mathrm{s}}\text{(4)}& {\displaystyle }{\mathrm{GPP}}_{max}=a{f}_{\mathrm{2}}{B}_{\mathrm{s}}\text{(5)}& {\displaystyle }\begin{array}{rl}{f}_{\mathrm{1}}\left(\mathrm{ds}\right)& =b{e}^{-c\mathrm{ds}}(\mathrm{ds}\ge d_{\mathrm{1}})\\ f_{\mathrm{1}}\left(\mathrm{ds}\right)& =g{e}^{h\mathrm{ds}}(\mathrm{ds}<d_{\mathrm{1}})\end{array}\text{(6)}& {\displaystyle }\begin{array}{rl}{f}_{\mathrm{2}}\left(\mathrm{ds}\right)& =i{e}^{-j\mathrm{ds}}(\mathrm{ds}\ge d_{\mathrm{2}})\\ f_{\mathrm{2}}\left(\mathrm{ds}\right)& =l\mathrm{ds}+m(\mathrm{ds}<d_{\mathrm{2}})\end{array}\end{array}$$

2.3 Brief introduction to the selected field site and information on the data for model evaluation

The field site (34^∘55.50^′ N, 110^∘42.59^′ E; altitude of 348 m) selected for this modeling case study is located at Dongcun Farm, near Yongji County, Shanxi Province, in northern China. The site is subject to a temperate continental monsoon climate, and it had an annual precipitation of 580 mm and a mean air temperature of 14.4^∘ from 1986 to 2010 (Cui et al., 2014). Cotton, winter wheat, and summer maize are the major crops grown at this farm and in the surrounding regions. Field experiments were performed in two adjacent fields (each of which was 100 m wide and 200 m long) for cotton and W–M between 2007 and 2010. The soil of the land cultivated with cotton and W–M was clay loam, with approximately 38% and 32 % clay (<0.002 mm), 57 % and 50 % silt (0.002–0.05 mm), 5 % and 18 % sand (0.05–2 mm), 10.0 and 11.3 g kg⁻¹ SOC, 1.1 and 1.1 g kg⁻¹ total nitrogen, a pH (in H₂O) of 8.0 and 8.7 at a depth of 0–10 cm, and bulk densities (0–6 cm depth) of 1.20 and 1.17 g cm⁻³, respectively (Liu et al., 2010, 2011, 2012). A sprinkler system was applied to both fields. For more detailed information on the field experiments and observed data, please refer to Cui et al. (2014), Liu et al. (2010, 2011, 2014, 2015), and Wang et al. (2013a, b).

The modified model was validated with observations from both fields. The daily meteorological data from 2004 to 2010 were obtained directly from Cui et al. (2014). The measured data were used directly for the minimum required soil properties. The input data on the field capacity and wilting point in WFPS were 0.65 and 0.2, respectively (Cui et al., 2014). The crop parameters for cotton were directly determined by the field measurements, which were 1900 kg C ha⁻¹ for potential grain yield (1.2 times the mean of the measured values), 0.41 and 25, 0.16 and 40, and 0.43 and 40 for the mass fractions and C∕N ratios of the grain, root, and leaf plus stem, respectively, and 3600 ^∘C for the TDD. Detailed management practices (Table S3) were obtained from Li et al. (2009) and Liu et al. (2014). In contrast to the locally conventional fertilizer application rate of 110–140 kg N ha⁻¹ yr⁻¹ for cotton, the fertilizer doses for the experimental cotton field in 2007 and 2008 were reduced to 66–75 kg N ha⁻¹ yr⁻¹ to avoid overgrowth of the leaves in place of seeds or lint. The data from the cotton field available for the model calibration and validation included the daily observed soil temperature (5 cm depth), topsoil moisture (0–6 cm), daily NEE from eddy covariance measurement from 2008 to 2009, sub-weekly observed CH₄ fluxes in the 2010 growing season (Liu et al., 2010, 2014, 2019; Wang et al., 2013b), daily N₂O and NO fluxes from 2007 to 2008, annually measured grain yields from 2008 to 2010, and annual ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching from 2008 to 2009. All the data used by Cui et al. (2014) in the validation of the original model were used to revalidate the modified model performance for the selected W–M fields with different field management practices. In addition, two observations of the cumulative NH₃ volatilization following urea topdressing to the winter wheat in the spring of 2008 (Tong et al., 2009) were also involved in the evaluation of the model performance (Yang et al., 2011). The information on all of the data used in the modified model calibration and validation is detailed in Table S1.

2.4 Scenario settings and simulations

For the three-crop system investigated, this study attempted to identify the BMP alternatives with the best rotation pattern between cotton and W–M and the optimized field management practices of this rotation pattern. For this purpose, six rotation pattern scenarios (hereinafter referred to as R₀, R₁, …, R₅) were set for a 6-year cycle that was widely applied by the local farmers (Liu et al., 2010, 2011, 2014). The subscript number of each rotation pattern represents the number of the consecutive years of cotton cultivation. For instance, R₀ denotes a 6-year monoculture of W–M, and R₂ represents a rotation pattern consisting of 2 years of continuous cotton rotated with 4 years of continuous W–M. Because the local farmers typically did not adopt a cotton monoculture for longer than 5 years, the longest cotton monoculture only lasted for 5 years (R₅). The transitions between cotton and W–M in the scenario rotations are detailed in Table S4.

As for the setting of the field management scenarios for the individual rotation patterns, four field management factors were considered, including (i) nitrogen fertilizer dose, (ii) water amount, (iii) method of irrigation, and (iv) depth of tillage. The values of these factors used for the baseline scenario were taken from the observations of the conventional management practices in the experimental region (Tables S3, S4, S5). The nitrogen doses of the baseline were 110 and 430 kg N ha⁻¹ yr⁻¹ for cotton and W–M, respectively. Over the last few decades, the fields in this region have been mostly flood-irrigated (Liu et al., 2010). Thus, flood-irrigation was chosen as the baseline method. The baseline timings and water amounts were established by referring to the 10–30 d cumulative precipitation prior to the individual irrigation events and the recorded timings and water amounts of the conventional management practices in the two adjacent fields. Thus, the irrigation frequencies and annual cumulative water amounts of the baseline varied from 1 to 3 times and 75 to 230 mm yr⁻¹ for cotton and 4 to 6 times and 290 to 510 mm yr⁻¹ for W–M (Table S5). The fraction of 100 % for residue incorporation and the conventional tillage to a depth of 20 cm for W–M and 30 cm for cotton were applied for the baseline practices. To screen the BMPs by fully taking the independent and interactive effects of rotation patterns and field management practices on the NIPs into account, a total of 6000 field management scenarios (each being a combination of the four management factors) for all the six rotation pattern scenarios (1000 for each) were randomly generated. The fertilizer doses and irrigation amounts were randomly selected, within their lower and upper bounds of continuous variations, using Latin hypercube sampling. The lower and upper bounds for nitrogen fertilizer doses (44–172 and 110–430 kg N ha⁻¹ yr⁻¹) and irrigation amount (40–100 mm per event) were set as 40 % and 100 % of the baseline, respectively. Only two irrigation methods, flooding (IF) – as the baseline – and sprinkling (IS), were included for the random sampling of this management factor. For cotton, the tillage depth was fixed at 30 cm. For W–M, four tillage depths (0 cm for no-till, 5 and 10 cm for reduced tillage, and 20 cm for the conventional practice) were included for random sampling of this factor. The BMPs for each rotation pattern scenario were first screened from 1000 field management scenarios. Then, the BMP alternatives were finally screened from the BMPs of the individual rotation pattern scenarios. These identified BMP alternatives would include the scenarios with overlapping NIP uncertainties, for which the random error at ±1 SD (instead of ±2 SD) for the total simulation error was adopted for the NIP of each scenario so as to achieve a relatively high screening precision. The decision variables and the NIP of baseline management scenarios for each rotation pattern scenario were used to particularly address the biogeochemical effects of rotation pattern.

An 18-year simulation was performed for each scenario. The annual averages for the simulated yields, decision variables, and NIPs were used to address the biogeochemical effects of rotation patterns or to screen the BMP alternatives. The simulations of all scenarios were driven by the meteorological data observed at Yuncheng station (approximately 60 km east of the experimental site) from 1996 to 2013 (http://data.cma.cn/data/cdcindex/cid/6d1b5efbdcbf9a58.html, last access: 10 October 2017). To stabilize the carbon and nitrogen dynamics and reduce the residual effects of the initial conditions (Palosuo et al., 2012; Zhang et al., 2015), a 12-year spin-up was performed for each scenario (i.e., a period of two consecutive 6-year rotation cycles) before the 18-year simulation. The spin-up for each scenario was driven by the same rotation pattern and field management practices as this scenario.

2.5 Method for identifying the best management practices

An objective method jointly relying on three constraints and NIPs was adopted in this study to screen the BMPs from given scenarios. These constraints included (i) stable or increased crop yields, (ii) an annual increase in the SOC stock by 4 ‰ or more, and (iii) reduced NEGE by 5 % or more. In the present study, the NEGE was determined by summing the emissions of CH₄ and N₂O as well as −ΔSOC, which was quantified as a CO₂ equivalent (CO₂eq) quantity. In the quantification of a NEGE, the global warming potentials at the 100-year time horizon, 34 for CH₄ and 298 for N₂O (IPCC, 2013), were used to convert the quantities of both gases into CO₂ equivalents. A NIP was expressed as a price-based proxy quantity in USD ha⁻¹ yr⁻¹ and used to evaluate the potential for a climatically and environmentally integrative impact exerted by a set of field management practices for multi-crop system. The NIP was determined as a linear function of the individual decision variables, following Eq. (7), wherein the multi-goal decision variables, NEGE, NH₃, NO, N₂O_ODM, and NL, represent the annual net ecosystem aggregate GHG emission (Mg CO₂eq ha⁻¹ yr⁻¹), NH₃ volatilization, NO emission, release of N₂O as an ozone layer depletion matter, and hydrological nitrogen loss (mainly by ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching), respectively (kg N ha⁻¹ yr⁻¹ for all the nitrogen-based variables). The coefficients k₁, k₂, k₃, k₄, and k₅ are mass-scaled price-based proxies for the NEGE, NH₃, NO, N₂O_ODM, and NL, with values of 7.00 USD Mg⁻¹ CO₂eq and 5.02, 25.78, 1.33, and 1.92 USD kg⁻¹ N, respectively (Cui et al., 2014). A lower NIP indicated a better set of management practices that can exert smaller negative impacts on the climate and environment. Accordingly, the BMP was identified as the scenario with the lowest NIP among the scenarios that could satisfy all three constraints.

$$\begin{array}{}\text{(7)}& \begin{array}{rl}\mathrm{NIP}& =k_{\mathrm{1}}\mathrm{NEGE}+k_{\mathrm{2}}{\mathrm{NH}}_{\mathrm{3}}+k_{\mathrm{3}}\mathrm{NO}+k_{\mathrm{4}}{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{ODM}}\\ & +k_{\mathrm{5}}\mathrm{NL}\end{array}\end{array}$$

2.6 Method for uncertainty quantification

The model simulation error (ε_s) of a NIP, a constraint variable (e.g., crop yield) or a decision variable involved in Eq. (7) represented the total simulation uncertainty. It was made of two components: the input-induced uncertainty (ε_input) due to the uncertainties of input items; and the model uncertainty (ε_model) due to insufficiencies in scientific structure or process parameters.

The ε_input of a simulated variable was a random error if the uncertainties of model input items were known as random errors. It was estimated using the Monte Carlo test with Latin hypercube sampling within the uncertainty ranges (95 % confidence interval – CI) of sensitive input items. In DNDC, the four basic soil properties (bulk density, pH, clay fraction, and SOC content) as input items were sensitive to the model outputs (e.g., Li, 2016). Accordingly, the uncertainties of these soil properties were regarded as the major sources of model output uncertainties, such as the constraint/decision variables, or NIP. The initial bulk density, pH, clay fraction, and SOC content ranged from 1.13 to 1.25 g cm⁻³, 8 to 8.7, 0.31 to 0.39, and 9 to 12 g kg⁻¹ (at the 95 % CI), respectively, which were adapted from the means and two times SD of spatially replicated observations in the two adjacent fields (Liu et al., 2014, 2019). A uniform distribution for each of these soil properties was assumed in the Monte Carlo test, in which the sampling and simulation were iterated until the mean of simulated NIPs for all iterations converged to a certain level within a tolerance of 1 %. The NIP uncertainty due to the model input uncertainties was presented as the SD of these iterated simulations.

An ε_s was systematic error reflecting a model simulation bias diverging from the truth. In this regard, the ε_s of a constraint/decision variable could be estimated using the slope of a zero-intercept univariate linear regression of simulations against observations (ZIR_s−o) or model relative biases (MRBs) resulting from model validation (Eq. 8). The MRBs were only used in case a significant ZIR_s−o was not obtained in model validation. To ensure a relatively high BMP screening precision, the random uncertainty of the ε_s of a NIP was presented as one times SD, instead of two times SD to represent the 95 % CI, as the BMP alternatives referred to the management scenarios with overlapping NIP uncertainties. For a constraint/decision variable, the mean or SD of ε_s in an absolute magnitude was estimated as the product of (i) an adjustment factor, (ii) the simulated variable quantity, and (iii) an error factor. The adjustment factor was obtained from model validation, which was estimated as the mean of the ratios of individual observations to simulations. The error factor for a variable with a significant ZIR_s−o was given as (mean-slope${}_{\mathrm{s}-\mathrm{o}}-\mathrm{1}$) ± SD-slope_s−o, wherein mean-slope_s−o and SD-slope_s−o denote the mean and SD of the ZIR_s−o slope, respectively. The item prior to and following the “±” was used to estimate the mean and SD of the ε_s. For a variable that failed to obtain a significant ZIR_s−o in model validation, the mean and SD of the error factors were given as the mean and SD of the MRBs. The mean of the ε_s for a NIP was estimated by simply summing the weighted absolute ε_s of individual decision variables. This was because the decision variables involved in the additive items of Eq. (7) were independent of each other. Meanwhile, the SD of the ε_s for a NIP was mathematically propagated from the SDs of the absolute ε_s of the decision variables.

In this study, the constraint variables included crop yield, −ΔSOC, and NEGE, although NEGE was also one of the decision variables. While there was no direct measurement of −ΔSOC and NEGE, their observation-oriented estimates were involved in the model validation, which provided the basis for the ε_s estimation of either variable. For the experimental fields with NEE observations, there was no significant input quantity of organic matter in manure or any other form while crop residues were fully incorporated into the soil. In these cases, each annual/seasonal −ΔSOC could be estimated as the sum of annual/seasonal NEE and yields, according to the mass conservation law, and used to represent the observation, with its uncertainty propagated from the random errors of the annual/seasonal yield and NEE measurements. The random error of this observation-oriented −ΔSOC that represented the annual/seasonal net CO₂ emission from the cropping system and those of the observed annual cumulative CH₄ and N₂O were propagated to estimate the observational error of the annual/seasonal NEGE. These observation-oriented estimates of −ΔSOC or NEGE were involved in model validation.

2.7 Statistics and analysis

Statistical criteria simultaneously used to evaluate the model validity included (i) the index of agreement (IA; Eq. 9), (ii) the NSI (Eq. 10; e.g., Moriasi et al., 2007; Nash and Sutcliffe, 1970), (iii) the determination coefficient (R²; Eq. 11) and the slope of a ZIR of observations against simulations (Jiang, 2010; Li et al., 2019), and (iv) the MRB (Eq. 8; e.g., Congreves et al., 2016; Willmott and Matsuura, 2005). In Eqs. (8)–(11), k and n (k=1, 2…, n) denote the kth pair and the total pair number of the values, respectively, $\stackrel{\mathrm{‾}}{o}$ represents the mean of the observations (o), and $\widehat{o}$ represents the predictions using the ZIR. The IA index fell between zero and one, with a value closer to one indicating a better simulation, and vice versa. An NSI value between zero and one indicated acceptable model performance. Better model performance was indicated by a slope and an R² value that were closer to one, and vice versa. An absolute MRB ($\left|\mathrm{MRB}\right|$) value smaller than the double coefficients of variation (CVs) of replicated observations implied a valid simulation (Dubache et al., 2019).

$$\begin{array}{}\text{(8)}& {\displaystyle }\mathrm{MRB}={\displaystyle \frac{s_k}{o_k}}-\mathrm{1}\text{(9)}& {\displaystyle }\mathrm{IA}=\mathrm{1}-{\displaystyle \frac{{\sum }_{k=\mathrm{1}}^n(s_k-o_k{)}^{\mathrm{2}}}{{\sum }_{k=\mathrm{1}}^n(\left|s_k-\stackrel{\mathrm{‾}}{o}\right|+\left|o_k-\stackrel{\mathrm{‾}}{o}\right|{)}^{\mathrm{2}}}}\text{(10)}& {\displaystyle }\mathrm{NSI}=\mathrm{1}-{\displaystyle \frac{{\sum }_{k=\mathrm{1}}^n(o_k-s_k{)}^{\mathrm{2}}}{{\sum }_{k=\mathrm{1}}^n(o_k-\stackrel{\mathrm{‾}}{o}{)}^{\mathrm{2}}}}\text{(11)}& {\displaystyle }{R}^{\mathrm{2}}=\mathrm{1}-{\displaystyle \frac{{\sum }_{k=\mathrm{1}}^n(o_k-{\widehat{o}}_k{)}^{\mathrm{2}}}{{\sum }_{k=\mathrm{1}}^n(o_k-\stackrel{\mathrm{‾}}{o}{)}^{\mathrm{2}}}}\end{array}$$

In this study, the statistical analysis and graphical comparison were performed using the SPSS Statistics client 19.0 (SPSS Inc., Chicago, USA) and Origin 8.0 (OriginLab, Northampton, MA, USA) software packages.

3 Results

3.1 Validation of daily simulations for both cropping systems

The seasonal dynamics and magnitudes of the soil temperature (5 cm) and topsoil moisture (0–6 cm) were predicted well by the model simulations (Fig. 1a, b). The sound model performance was indicated by the IA, NSI, and ZIR slope and R² values of 0.98 and 0.83, 0.93 and 0.15, 0.93 and 0.83, and 0.95 (n=677, p<0.001) and 0.42 (n=432, p<0.001) for the temperature and moisture, respectively.

Figure 1 Observed and simulated daily mean soil temperature (5 cm), soil moisture (0–6 cm), daily net ecosystem exchanges of carbon dioxide (NEE) in cotton field and winter wheat–summer maize fields, and daily fluxes of methane (CH₄), nitrous oxide (N₂O), and nitric oxide (NO) from the cotton field. The solid and dashed arrows indicate the dates of fertilization and irrigation, respectively. The measurement errors were not shown in panels (a)–(e) for figure clarity. The vertical bar for each observation in panels (f)–(h) indicates ±2 SD to represent the uncertainty at the 95 % confidence interval. The legend in (c) applies to all panels.

Compared with the original model, the modified model showed much better performance in simulating daily NEE fluxes from both cropping systems. It simulated the abrupt NEE fluxes on cloudy or rainy days in the growing season particularly well (Fig. 1c–e). In comparison with the original model for daily NEE simulations, the modified model enhanced the IA, NSI, and ZIR slope and R² values from 0.74 to 0.81, 0.32 to 0.60, 0.60 to 0.96, and 0.27 to 0.60 (n=261, p<0.001), respectively, for the cotton field, and from 0.75 to 0.80, 0.30 to 0.51, 0.69 to 0.80, and 0.45 to 0.55 (n=311, p<0.001), respectively, for the W–M field. For the CH₄ uptake, the observations and simulations showed similar seasonal variations (Fig. 1f), with IA, NSI, and ZIR slope and R² values of 0.68, <0, 0.70, and 0.08 (n=69, p=0.018), respectively.

The simulated seasonal patterns and peak emissions of N₂O and NO generally agreed with the observations (exemplified by Fig. 1g, h for the cotton field). Compared with the original model for daily N₂O flux simulations, the modified model performed comparably well for the cotton field, with a NSI of around −0.45, and ZIR slope ∼0.49 and R² values of ∼0.39 (n=592, p<0.001), although it showed better performance for the W–M fields, with IA, NSI, and ZIR slope and R² values of 0.69 versus 0.52, 0.03 versus –0.26, 0.62 versus 0.46, and 0.16 (n=976, p<0.001) versus “not available”, respectively. Relative to original model, the modified model showed improved simulations for the daily NO fluxes from the cotton field, with increased IA, NSI, and ZIR slope and R² values from 0.62 to 0.78, −1.03 to −0.04, 0.37 to 0.54, and 0.09 to 0.39 (n=333, p<0.001), respectively, while it performed comparably well for those from the W–M fields, with IA, NSI, and ZIR slope and R² values of ∼0.82, ∼0.32, ∼0.74, and 0.35–0.40 (n=967, p<0.001), respectively.

In comparison with the observed daily NH₃ fluxes (measured using micrometeorological technique) following a single fertilization event in the maize season in 2008, the modified model simulations showed IA, NSI, and ZIR slope and R² values of 0.87, 0.12, 0.68, and 0.53 (n=11, p=0.07), respectively. However, the model failed to capture the immediate responses of daily NH₃ fluxes to urea addition in the wheat fields, which were measured using a quasi-dynamic chamber method.

3.2 Validation of simulated variables in annual/seasonal cumulative quantities

For the cotton yields (seeds plus lint) over the three consecutive experimental years (2008–2010), the simulations were consistent with the observations in terms of the smaller $\left|\mathrm{MRBs}\right|$ (0.4 %–24 %) than the double CVs (39 %–56 %) of the spatially replicated measurements. For all of the experimental treatments of the cotton, wheat, and maize, the simulated yields of both the modified and original model strongly agreed with observations (Fig. 2a), with an IA of 0.93–0.95, a NSI of 0.75, and ZIR slope and R² values of 0.96–1.00 and 0.75–0.78, respectively (n=35, p<0.001). This validation resulted in an adjustment factor of 0.96 and smaller error factors of 3.0±1.6 % for crop yields simulated by the modified model.

Figure 2 Comparison between observations and simulations of crop yields, annual/seasonal cumulative NEE and NEGE, annual/seasonal ΔSOC, annual cumulative fluxes of methane (CH₄) uptake, nitrous oxide (N₂O) and nitric oxide (NO), and cumulative fluxes of ammonia (NH₃). Yield refers to the seed yield of cotton (open circle) and the grain yield of winter wheat (solid circle) and summer maize (solid diamond). NEE refers to the net ecosystem exchanges of carbon dioxide. NEGE refers to the net ecosystem aggregate greenhouse gas emission. ΔSOC refers to the change in soil organic carbon stock. NEE, NEGE, and ΔSOC are annual for cotton and seasonal for wheat and maize. The observed ΔSOC was given as the opposite of NEE plus yield in carbon mass quantity for the cropping system with the incorporation of full residues, whereas each ΔSOC simulation was the sum of simulated changes in carbon stocks of soil humus, microbial biomass, and dissolvable organic compounds. Simulations resulted from the modified model. The slope errors of the zero-intercept linear regressions are ±2 SD to represent the 95 % confidence interval. Vertical bars indicate the standard deviation of three or four spatial replicates, with the exception of NEE. The errors of NEE were adapted from the coefficient of variation (on average 25 %) reported by Wang et al. (2013b). DM refers to dry matter. CO₂eq refers to the carbon dioxide equivalent. The 100-year global warming potentials of 34 for CH₄ and 298 for N₂O (IPCC, 2013) were used to quantify NEGE in CO₂eq quantity.

For the annual cumulative NEE of the cotton field during the two consecutive year-round periods and the seasonal cumulative NEE in two wheat seasons and one maize season, the simulations of both model versions showed comparably significant agreement with the observations (Fig. 2b), with an IA of 0.99–1.00, a NSI of 0.95–1.00, and ZIR slope and R² values of 0.92–1.02 and 0.97–0.99, respectively (n=5, p≤0.000–0.002). The modified model simulations showed $\left|\mathrm{MRBs}\right|$ of 6 %–16 %, which were much less than the reported CV (25 %) of the eddy covariance observations.

Compared with the annual/seasonal NEGE quantities derived from the observations of crop yields, annual/seasonal cumulative NEE, and fluxes of CH₄ and N₂O, the simulations implicated good performance of the modified model (Fig. 2c), with an IA of 0.96, a NSI of 0.77, and ZIR slope and R² values of 0.73 and 0.91 (n=5, p=0.013). Although the simulations showed an average overestimation of ∼25 %, their $\left|\mathrm{MRBs}\right|$ were only 17 %–72 % (33 % on average) of the observation-oriented CVs (27 %–170 %, with a mean of 91 %), implicating a statistically meaningful good performance of the modified model. This validation resulted in an adjustment factor of 0.73 and error factors of 25±19 % for the annual/seasonal NEGE simulated by the modified model.

In comparison with the annual/seasonal ΔSOC quantities estimated from the observed crop yields and annual/seasonal cumulative NEE, the simulations by the modified model (Fig. 2d) showed an IA of 0.96, a NSI of 0.75, and ZIR slope and R² values of 0.71 and 0.92 (n=5, p=0.011). The ZIR slope indicated an average overestimation of the model by ∼30 %. Nevertheless, the $\left|\mathrm{MRBs}\right|$ of the individual simulations were only 4 %–79 % (30 % on average) of the observation-oriented CVs (30 %–210 %, with a mean of 97 %), still indicating a statistically meaningful consistency. This validation resulted in an adjustment factor of 0.71 and error factors of 30±19 % for the annual/seasonal ΔSOC simulated by the modified model.

The model simulations of the annual cumulative CH₄ uptake in 2009 and 2010 showed significant agreement (Fig. 2e), with an IA of 0.98, a NSI of 0.91, and ZIR slope and R² values of 1.00 and 0.91 (n=7, p<0.001). The $\left|\mathrm{MRBs}\right|$ were only 5 %–56 % (25 % on average) of the double CVs (10 %–24 %, with a mean of 17 %) for the spatially replicated observations. This validation resulted in an adjustment factor of 1.00 and error factors of $-\mathrm{0.2}\pm \mathrm{1.7}$ % for the modified model simulations of cumulative CH₄ uptake.

The modified model simulations of the annual cumulative N₂O emissions from all of the field experimental treatments of the cotton and W–M fields were comparable with the observations (Fig. 2f), showing an IA of 0.94, a NSI of 0.72, and ZIR slope and R² values of 0.90 and 0.83 (n=12, p<0.001). The $\left|\mathrm{MRBs}\right|$ of the individual simulations were only 6 %–93 % (36 % on average) of the double CVs (23 %–64 %, with a mean of 47 %) for the spatially replicated observations. For the annual cumulative N₂O emissions simulated by the modified model, this validation resulted in an adjustment factor of 0.90 and error factors of 8.7±4.5 %.

Compared with annual cumulative N₂O emissions, slightly better consistency with observations was obtained for the modified model simulations of the annual cumulative NO emissions from the cotton and W–M fields under different experimental conditions (Fig. 2f, g). The NO simulation showed an IA of 0.97, a NSI of 0.85, and ZIR slope and R² values of 0.90 and 0.94 (n=11, p<0.001). The $\left|\mathrm{MRBs}\right|$ of the individual simulations were only 2 %–52 % (23 % on average) of the double CVs (30 %–99 %, with a mean of 66 %) for the spatially replicated observations. This validation provided an adjustment factor of 0.90 and error factors of 10.1±3.2 % for the cumulative NO emissions simulated by the modified model.

The simulations of cumulative NH₃ volatilization during the 11 d following the three urea application events, with one in the maize field in summer and two in the winter wheat fields (with flood-irrigation and sprinkling, respectively) in spring (Fig. 2h), showed an IA of 0.97, a NSI of 0.86, and ZIR slope and R² values of 1.00 and 0.86 (n=3, p=0.246). The simulations resulted in smaller MRBs of −8.8 % to 3.8 % (with a mean of −0.4 %), despite the model's failure to capture the quick responses of daily NH₃ fluxes to the urea topdressing events. This validation resulted in an adjustment factor of 1.00 and error factors of $-\mathrm{0.4}\pm \mathrm{7.3}$ % for the modified model simulations of cumulative NH₃ volatilization following nitrogen applications.

The modified model simulations of the cumulative ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching from the cotton field in 2 consecutive years agreed with the observations, in terms of the smaller MRBs (−32 % to −27 %) than the two times CVs (109 %–115 %) for the spatially replicated observations. These MRBs represented the model underestimations by 3–4 and 13–21 kg N ha⁻¹ yr⁻¹ for the annual ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching rates in the cotton and W–M fields subject to the currently applied field management practices, respectively. This validation derived an adjustment factor of 1.42 and error factors of $-\mathrm{29}\pm \mathrm{4}$ % for the modified model simulations.

The above results suggested that the modified DNDC95 model was especially applicable at this field site for investigating the biogeochemical effects of different rotation patterns between the cotton and W–M and those exerted by different management practices, and thus was capable of BMP identification.

3.3 Biogeochemical effects of different cotton and wheat–maize rotation patterns

Figure 3 illustrated the dynamics of the crop yields and each decision variable resulting from the consecutive simulations over 18 years for all of the rotation patterns subject to the field management practices of the baseline scenario. Figure 4 showed the relationship between the annual average of each decision variable and the number of consecutive years of cotton monoculture within the rotation patterns.

Figure 3 Simulated cumulative crop yields, changes in soil organic carbon (ΔSOC), methane (CH₄), nitrous oxide (N₂O) releases, net ecosystem aggregate greenhouse gas emission (NEGE), ammonia (NH₃) volatilization, nitric oxide (NO) emission, and nitrate leaching (NL) of individual rotation patterns (with a 6-year rotation cycle) over an 18-year period. R₀, R₁, …, R₅ represent rotation patterns with the cotton cultivated consecutively for 0, 1, …, 5 year(s), respectively, within each 6-year rotation cycle. The legend in (e) applies to all panels. The simulations resulted from the modified model driven by the currently applied field management practices (i.e., the baseline field management scenario) and the observed means of input soil properties.

The average grain yields for the cotton, wheat, and maize were not significantly different among the various rotation patterns, with averages of 3.5, 4.8, and 6.7 kg dry matter ha⁻¹ for cotton, wheat, and maize, respectively (Fig. 3a–c).

For the dynamic changes in the annual SOC stocks, the values were generally positive for W–M but negative for cotton, except for the first year after the transition to this fiber crop. As indicated by Fig. 3d, the simulated SOC contents over the 18-year period increased for R₀, R₁, R₂, and R₃ but decreased for R₄ and R₅. The annual average −ΔSOC increased significantly (p<0.001) with an increase in the consecutive years of cotton monoculture from 0 to 5 within the 6-year rotation cycle (Fig. 4a). The rotation patterns with the baseline management showed small variations in the CH₄ uptake (Fig. 3e), with the uptake rates ranging from 1.6 to 2.1 kg C ha⁻¹ yr⁻¹. However, the annual averages for the CH₄ uptake increased significantly (p<0.001) with the increased consecutive years of cotton monoculture (Fig. 4b). For N₂O, the annual emissions showed large inter-annual variations (Fig. 3f), with CVs of 26 %–48 %. In addition, the average emissions of this gas decreased significantly from 4.6 to 2.6 kg N ha⁻¹ yr⁻¹ (Fig. 4c) after increasing consecutive years of cotton monoculture (p<0.001). As a result, the NEGE was significantly promoted (p=0.002) (Figs. 3g, 4d).

Figure 4 Simulated effects of various rotation patterns between a cotton and winter wheat–summer maize cropping system with a 6-year cycle on decision variables and the negative impact potential (NIP). The subscript numbers in R₀, R₁, …, R₅ refer to the number of consecutive years of cotton cultivation. The y axis units are megagrams of carbon per hectare per year (Mg C ha⁻¹ yr⁻¹) for the opposite of the mean annual increase in soil organic carbon stock (−ΔSOC), kilograms of carbon per hectare per year (kg C ha⁻¹ yr⁻¹) for methane (CH₄) emissions, kilograms of nitrogen per hectare per year ( kg N ha⁻¹ yr⁻¹) for fluxes of nitrous oxide (N₂O), ammonia (NH₃), nitrous oxide (NO), and nitrate leaching (NL), megagrams of carbon dioxide equivalent per hectare per year (Mg CO₂eq ha⁻¹ yr⁻¹) for net ecosystem aggregate greenhouse gas emission (NEGE), and US dollars per hectare per year (USD ha⁻¹ yr⁻¹) for NIP. The CO₂eq was based on the 100-year global warming potentials, i.e., 34 for CH₄ and 298 for N₂O (IPCC, 2013). The NIP was calculated using Eq. (7), which is presented in the text. The vertical bar within the open circle of each datum point indicates the absolute uncertainty (1 SD) introduced by the input uncertainties of key soil properties. Each unfilled column indicates the absolute total uncertainty of the simulation, with its vertical bar representing its random uncertainty (1 SD).

Regarding the gaseous air pollutants NH₃ and NO, the simulated emissions ranged from 17 to 103 and 0.5 to 3.3 kg N ha⁻¹ yr⁻¹, respectively (Fig. 3h, i). Figure 4e and f showed that the average annual emissions of both gases were significantly reduced after increasing the consecutive years of cotton monoculture (p≤0.001). The annual ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching of the different rotation patterns displayed significant inter-annual variations (Fig. 3j), with CVs of 41 %–69 %. Thus, the annual averages for ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching changed insignificantly in response to the consecutive years of cotton monoculture (p=0.056; Fig. 4g).

The NIP varied significantly among the various rotation patterns (p<0.001), declining from 610 to 324 USD ha⁻¹ yr⁻¹ with increased consecutive years of cotton monoculture (Fig. 4h). For the three constraints, the crop yields showed no obvious differences among the various rotation patterns. Both R₀ and R₅ represented the typical rotation patterns in the region. The simulations for the former indicated the greatest increase in SOC and the lowest NEGE but the highest NIP, whereas those for the latter showed the greatest SOC loss and the largest NEGE but the lowest NIP (Fig. 4a, d, h). These patterns indicated that neither typical rotation pattern is sustainable.

3.4 Identification of best management practices

Out of the 6000 field management scenarios, three under the R₃ were finally identified as the BMP alternatives, which simultaneously satisfied the given constraints while yielding the comparably lowest NIPs (332–335 USD ha⁻¹ yr⁻¹) within the overlapping uncertainty ranges with ε_s of $-\mathrm{22}\pm \mathrm{16}$ USD ha⁻¹ yr⁻¹ (Table 1). These BMP alternatives for the three-crop system recommended the following combination of field management practices: (i) the currently applied 6-year rotation cycle with 3-year cotton monoculture rotated with 3-year W–M; (ii) full incorporation of crop residues at harvest; (iii) the presently adopted crop cultivars and timing of sowing, fertilization (date, depth, and splits), irrigation (date and times), and harvest; (iv) the use of urea alone and 18 % lower rates (90 and 353 kg N ha⁻¹ yr⁻¹ for the cotton and W–M, respectively) than the conventional nitrogen fertilization; (v) sprinkling or flood-irrigation using ∼23 % less water (∼77 mm per event) than conventional flood-irrigation; and (vi) conventional plow tillage (30 cm depth) following final cotton harvest but reduced tillage (rotary 5 cm depth) for the W–M. In comparison with the simulations driven by the baseline scenario (R₃ as the currently applied rotation pattern and its field management practices), the identified BMP alternatives could sustain the baseline crop yields while simultaneously enlarging the SOC stock by ≥4 ‰ and mitigating the NEGE, NH₃ volatilization, NO emission, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching by ∼7 %, ∼25 %, ∼2 %, and ∼43 %, respectively, despite a slight increase (by ∼5 %) in the N₂O emission (Table 1).

Table 1 Simulated constraint and decision variables and negative impact potentials (NIPs) for the baseline (the conventionally applied practices) and the alternatives of the best management practices.

BAS: the baseline. BMP: different best management alternatives denoted by subscript numbers. N: nitrogen fertilizer dose (kg N ha⁻¹ yr⁻¹) of cotton/wheat–maize (W–M). IrA: irrigation amount (mm per event). IM: irrigation method. IF: flood-irrigation.
IS: sprinkling irrigation. T: tillage depth (cm). R: rotation pattern. R₃: rotation pattern with 3 consecutive years of cotton rotated with 3 continuous years of W–M. Yield: seed (cotton) or grain (W–M) yield (Mg ha⁻¹ in dry matter). Sim: annual quantity simulation.
ε_s: the absolute total simulation error (i.e., the systematic error) of the annual quantity simulation, with its error (1 SD) representing the random uncertainty magnitude. ε_input: the random model simulation error (1 SD) due to input uncertainties of the key soil
properties including clay fraction, bulk density, pH, and soil organic carbon content. ΔSOC: annual change in soil organic carbon stock at a depth of 0–50 cm (Mg C ha⁻¹ yr⁻¹). NEGE: the net ecosystem aggregate greenhouse gas emissions in carbon dioxide
equivalent (Mg CO₂ eq ha⁻¹ yr⁻¹). CH₄: methane emissions (kg C ha⁻¹ yr⁻¹). N₂O, NH₃, and NO: the emissions of nitrous oxide, ammonia, and nitric oxide, respectively (kg N ha⁻¹ yr⁻¹). NL: nitrate leaching (kg N ha⁻¹ yr⁻¹). The CO₂eq was based on
the 100-year global warming potentials of 34 for CH₄ and 298 for N₂O (IPCC, 2013). NIP: negative impact potential (USD ha⁻¹ yr⁻¹). The simulations are the averages of 18 consecutive years.

4 Discussion

4.1 Model performance

The DNDC model has been widely applied in agricultural systems around the world. The version modified in this study showed good performance in simulating the soil environmental factors (soil temperature and moisture), crop yields, NEE, NH₃ volatilization, CH₄ uptake, emissions of N₂O and NO, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching for the investigated lands cultivated with cotton and W–M under different field management treatments. The satisfactory validations of both crop systems, especially for the constraint and decision variables at the annual scale, suggested that the modified DNDC95 could be applied to quantify the constraint and decision variables to determine the NIP for the cotton and W–M rotation system under various management practices.

The well-simulated soil environmental factors and crop yields provided a solid basis for further simulating the constraint and decision variables under any field management condition of the three-crop rotation system. This is because the soil environmental factors are the key factors regulating the biogeochemical processes, and crop yields are indicators of essential processes in plant nitrogen uptake (Chirinda et al., 2011; Kröbel et al., 2010). For the simulations of the N₂O and NO emissions, discrepancies in daily emissions generally occur in the simulations of DNDC or other current biogeochemical models due to the interactions among soil environmental factors and complex carbon- and nitrogen-related processes (e.g., Bell et al., 2012; Chirinda et al., 2011; Cui et al., 2014; Lehuger et al., 2011), which may occasionally result in significant time lags between the observations and simulations (e.g., Zhang et al., 2015). For the cotton in this study, the significant underestimation of daily NO fluxes in spring was solved to some extent by modifying the model version used by Cui et al. (2014). However, this improvement did not significantly affect the annual cumulative emissions, which were not primarily contributed by the spring fluxes (Liu et al., 2015). In fact, occasional time lags of one to a few days for measured/simulated daily fluxes seldom lead to a significant modification for the seasonal/annual cumulative emissions of a nitrogenous gas. This is attributed to the control of the mass conservation law and the canceling effect of negative and positive daily errors. The modified algorithm improved the simulations of daily NEE fluxes, thereby providing a solid basis for yielding reliable annual/seasonal cumulative NEE quantities. According to the mass conservation law, the annual cumulative NEE can be included as one of the two additive items to estimate the annual ΔSOC of an annual crop system with retention/incorporation of full residues but without significant input/output of organic matter other than product removal at harvest. This approach may be used as an alternative algorithm in the modified model to simulate the annual ΔSOC of such a cropping system. In the present case study, the annual ΔSOC simulated by the modified model using this alternative approach was consistent with those simulated by the algorithm that quantifies the annual ΔSOC by summing the annual carbon pool changes in the humus, microbial biomass, and dissolvable organic compounds (Cui et al., 2014). The consistency was indicated by the $\left|\mathrm{MRBs}\right|$ of 20±50 % versus 37±117 % (95 % CI) in comparison with the three observation-oriented estimates of annual ΔSOC (two for cotton and one for the W–M). Due to the marginally small sample size (n=3), this preliminary result still requires confirmation through further study. The simulated NH₃ volatilization from the cotton field accounted for 18 %–24 % of the applied fertilizer nitrogen during the two year-round periods involved in the model validation. These simulated nitrogen loss rates via NH₃ volatilization were comparable with the reported field measurements of 10 %–23 % (Li et al., 2016).

The model validation in this study suggested that the satisfactory simulations of constraint and decision variables at the annual scale could provide a solid basis for BMP identification. Because of the limited annual observations of NH₃ volatilization, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching, and ΔSOC estimated from annually measured NEE, the insufficient validation still resulted in large uncertainties in the simulations of these three variables. Therefore, future studies are still required for further validation of the model performance using comprehensive observations covering these variables as well as the others, thus reducing the simulation errors of the constraint and decision variables so as to improve the screening precision of BMP alternatives.

4.2 Biogeochemical effects of the rotation pattern and other management practices

The scenario analysis relying on model simulations in this study showed that environmental nitrogen contamination could be reduced while (i) sustaining crop yields to protect food security, (ii) achieving the 4 ‰ goal in soil carbon sequestration, and (iii) decreasing the net ecosystem aggregate GHG emissions to mitigate climate change. The reductions in environmental nitrogen contamination could be attributed to the better synchronization of crop nitrogen requirements and soil nitrogen availability by optimizing field management practices.

For cotton, a period of 5 consecutive years is usually applied as the longest cotton monoculture to stabilize its yields. Within this period, balanced elemental nutrients are applied, and thus the negative effect of monoculture on cotton yields can be offset in practice (Han, 2010). In addition, the DNDC model assumes balanced nutrient supplies for any crops as well as optimum phytosanitary conditions; thus, the negative effects of monoculture are not taken into account (e.g., Li, 2016).

The simulated positive annual ΔSOC for the W–M cropping system were mainly attributed to the incorporation of the full aboveground residues (at rates of 5.1–7.0 Mg C ha⁻¹ yr⁻¹), which in turn favored carbon sequestration (Han et al., 2016). On the contrary, the simulations of annual ΔSOC for the cotton cropping system were negative. The SOC stock decreases resulted from (i) the more notable CO₂ emissions over the longer fallow season and (ii) the lower rates of fully incorporated residues (at rates of 2.5–3.1 Mg C ha⁻¹ yr⁻¹) than those of W–M system (Liu et al., 2019). As a remarkable carbon sink, the W–M cropping system with the incorporation of the full crop residues could even completely compensate for the SOC lost during the first cotton-planting year following the transition. Thus, the simulated annual ΔSOC was generally positive during the first cotton cultivation year of a three-crop rotation cycle. As a result, the R₀ (i.e., pure W–M continuously within each 6-year period) acted as a net GHG sink; this was due to the fact that the positive ΔSOC could exceed the N₂O emission during the W–M cultivation, whereas all of the three-crop systems subject to R₁–R₅ rotation patterns would function as net GHG sources. The higher nitrogen application rate for the W–M system than for the cotton system resulted in more reactive nitrogen remaining in the soil (Chen et al., 2014; Ju et al., 2009), thereby stimulating higher N₂O and nitrogenous air pollutant emissions in the trials with fewer cotton cultivation years. Therefore, the rotation patterns of cotton and W–M can be optimized to realize sustainable intensification in terms of sustaining crop yields at a relatively high level, maximizing SOC increase, and minimizing negative impacts on the climate and environment.

Northern China, as the most important agricultural region in the country, experienced an increase in crop yields by a factor of 2.8 from 1980 to 2008. During this period, the application of mineral fertilizers increased by a factor of 5.1. The rapid increase in fertilizer use has resulted in excessive nitrogen remaining in the soil, posing potential risks for the environment (Chen et al., 2011; Zhang et al., 2017). To solve this problem, a reduction in fertilizer application has been proposed in several previous studies (e.g., Chen et al., 2011, 2014; Liu et al., 2012). The results of the scenario analysis in this study indicated that further reducing the farmer-optimized nitrogen doses by 18 % could still sustain the crop yields while greatly decreasing the release of nitrogenous pollutants.

In addition to fertilization, over-irrigation has also been ubiquitous in northern China for a long time, and is threatening the water security of this region due to the sharply declining groundwater table and water pollution (Gao et al., 2015; Ju et al., 2009). For this reason, only management options that can reduce the amount of irrigation should be recommended due to the severe shortage of water resources in the region. In addition, adopting sprinkling irrigation instead of flood irrigation for an equal amount of water showed positive effects on the crop yields, indicating improved irrigation efficiency (Zhang et al., 2017). This result indicated that increasing the water-use efficiency via the application of alternative irrigation techniques coupled with reduced nitrogen doses could be a pathway to sustain crop yields.

Reduced tillage practices have been promoted in China in the recent years. To facilitate the decomposition of the woody cotton residues and avoid outbreaks of diseases and pests induced by continuous implementation of reduced tillage or no-till, the tillage practices were only adjusted in the W–M fields, whereas the currently applied deep tillage practices were maintained for the cotton when setting the tillage scenarios. The simulations showed that the reduced tillage or no-till practices could sustain the crop yields while also reducing NH₃ volatilization and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching; these results were consistent with the reports from other experimental studies (e.g., Zhao et al., 2016).

As shown above, the appropriate combinations of a rotation pattern and field management practices can satisfy the three given constraints while resulting in the lowest NIPs with overlapping uncertainties. However, direct observations in field experiments usually with very limited management treatments are far less sufficient for screening these appropriate combination alternatives. Nevertheless, identifying the appropriate combination alternatives is one of the purposes of biogeochemical models, such as DNDC. In principle, a biogeochemical model that is validated with limited observations from field experiments, such as the DNDC95 modified and used in this study, could be capable of fulfilling this task.

4.3 Evaluation of the best management practice

The scenario analysis in this study was effective for screening the BMP alternatives. The identified BMP alternatives could sustain the crop yields of the three-crop rotation system, increase the SOC stock annually at 4 ‰ or more, mitigate the NEGE, and reduce the NH₃ and NO emissions as well as ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching due to the enhanced resource use efficiency in response to the reduced nitrogen-fertilizer doses, irrigation amounts, and tillage depth for the W–M system. Hence, the BMP alternatives could result in significantly reduced NIPs even compared to the currently applied field management practices that have been optimized by the local farmers. However, the BMP alternatives identified were based on the constraint and decision variables only validated against the observations at the single field site involved in this case study. In this regard, confirmation of these BMP alternatives at other sites in this region is still required in the future studies.

A biogeochemical model as an ideal tool for identifying the BMPs is reflected by near-zero ε_s for any constraint/decision variable or NIP. A small sample size of the observations used for validation of any constraint/decision variable would result in largely positive or negative ε_s (including large over- or under-estimations) for model simulations of the variable, likely account for the large ε_s of a NIP, and thus lead to a lower precision in screening the BMP alternatives. Therefore, the applicability of the approach proposed in this study for identifying the BMPs is highly dependent upon the validations using observations with appropriate sample sizes for individual constraint and decision variables. In this study, the ε_s and ε_input for the simulated variables and NIPs of management scenarios were quantified. For the NIPs of the identified BMP alternatives, for instance, the ε_s and ε_input at relative magnitudes were 6.5±4.9 % and ±3.3 %, respectively, which were similar to those (9.1±5.0 % and ±3.1 %, respectively) of the baseline scenario. According to these errors, the uncertainty ranges of the NIPs for the three alternatives almost fully overlapped with each other, and they were all beyond the uncertainty range of the NIP for the baseline scenario. This implied that the approach proposed in this study could be applicable for identifying the BMPs of the three-crop rotation system. Nevertheless, the ε_s uncertainty range of one times SD still fully diverged negatively from zero, due to the marginally small sample sizes of available ΔSOC, NH₃ volatilization, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching observations that led to insufficient validations for these variables. Especially, the model underestimations of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching (with an adjustment factor of 1.42 and error factors of $-\mathrm{29}\pm \mathrm{4}$ %) overwhelmingly dominated the diverged ε_s of the NIPs, which were comparable to or larger than the ε_input values. Relying on the few field observations, one was still not able to judge whether there were insufficiencies in the scientific structures or inappropriate parameters in the model to dominate the large ε_s for these variables. Therefore, multiple comprehensive field observations with appropriate sample sizes to fully cover all the relevant variables are necessary for use in an advanced biogeochemical model with multiple functions in order to address the best management issue of a multi-crop rotation system to achieve multiple benefits.

The DNDC model was established by following the mass conservation law. In other words, this model can accurately reflect the mass balance of the carbon or nitrogen budgets for the simulated soil layer (0–50 cm depth). This principle implies that only one nitrogen budget item could be omitted for model validation. This item is usually soil nitrogen loss via the production of dinitrogen gas (N₂), mainly by denitrification, which is very difficult to measure in situ (e.g., Wang et al., 2013; Zhang et al., 2019). For both crop cropping systems, however, the nitrogen lost via this pathway could almost be fully inhibited in the topsoil, wherein the soil moisture contents were often lower than 60 % WFPS (Linn and Doran, 1984; Liu et al., 2011, 2014). For instance, the N₂ emission likely accounts for approximately 1.6 % of the urea applied in a winter wheat season (Zhang et al., 2019), which is at a negligibly low level for the nitrogen balance.

Regarding the identification of the BMPs, the approach proposed and applied in this case study only includes the biogeochemical effects of management on the constraint/decision variables. This approach currently excludes other factors, such as those related to the costs of the management practices, thereby likely resulting in uncertainties in the screened BMP alternatives. Despite some deficiencies, this approach can be easily and automatically implemented as long as the simulations for all constraint and decision variables can be validated using comprehensive observations, implicating its potential applicability for more comprehensive situations. Adding the missing factors is one of the future research tasks to further improve this approach.

5 Conclusions

To address the challenging issue of optimizing multi-crop system management to simultaneously achieve multiple benefits, a biogeochemical model-based approach for identifying the best management practices (BMPs) was proposed and tested in this site-scale case study. A widely distributed three-crop system in northern China, which grew cotton in rotation with winter wheat and summer maize (W–M), was investigated. The BMPs referred to the management alternatives with the lowest negative impact potentials (NIPs), falling within overlapping uncertainty ranges, among the scenarios satisfying a set of constraints. The NIP of a scenario was defined as the linear function of five decision variables, including the net ecosystem aggregate greenhouse gas emission (NEGE), ammonia (NH₃) volatilization, nitric oxide (NO) release, emissions of nitrous oxide (N₂O) as an ozone layer depletion matter, and nitrate leaching. This study used three variables, i.e., crop yield, annual change in SOC stock (ΔSOC), and NEGE, to specify the applied constraints that were stable/increased crop yields, annual ΔSOC by 4 ‰ or more, and reduced annual NEGE by at least 5 % in comparison with those of the baseline scenario (as the currently applied practices in this study). The constraint and decision variables to determine the NIP of each scenario were provided by a simulation using the DeNitrification-DeComposition version 95 model (DNDC95), which was modified in this study. Due to the unsatisfactory performance of the model in daily simulations of NO emission and net ecosystem exchange of carbon dioxide (NEE), the model was modified to include a new parameterization of soil moisture effects on the NO production during nitrification and replacement of the original calculation approach for NEE with an algorithm based on gross primary production. For the variables with available measurements in two adjacent fields at the selected study site, the modified model showed statistically meaningful consistency between simulations and observations. Using the systematic errors obtained from the model validation to determine the simulation uncertainties of the concerned variables for each scenario and that of its NIP, the modified model simulations driven by 6000 management scenarios automatically identified three BMP alternatives. These BMP alternatives follow the current adopted rotation pattern (3 consecutive years of cotton rotated with 3 continuous years of W–M), with 18 % less fertilizer nitrogen, ∼23 % less irrigation via sprinkling or flood-irrigation, and reduced depth of tillage for the W–M in comparison with the current applied farmer-optimized management practices. This case study demonstrated the practicability of the model-based approach and outlined its potential applicability for optimizing the field management of a multi-crop system to simultaneously achieve multiple United Nations Sustainable Development Goals. It also emphasized the need to make comprehensive observations that fully cover the constraint and decision variables, other related factors, as well as all of the crops and management practices in question to facilitate effective BMP screening through virtual experiments using a biogeochemical model, such as DNDC. In future study to identify the BMPs, specifically for the three-crop rotation system at the regional scale, it is still necessary to use a 6-year model validation that includes a rotation of all three commodity crops as well as all of the studied management practices in question.