A Sensitivity Analysis of the SPACSYS Model

: A sensitivity analysis is critical for determining the relative importance of model parameters to their inﬂuence on the simulated outputs from a process-based model. In this study, a sensitivity analysis for the SPACSYS model, ﬁrst published in Ecological Modelling (Wu, et al., 2007), was conducted with respect to changes in 61 input parameters and their inﬂuence on 27 output variables. Parameter sensitivity was conducted in a ‘one at a time’ manner and objectively assessed through a single statistical diagnostic (normalized root mean square deviation) which ranked parameters according to their inﬂuence of each output variable in turn. A winter wheat ﬁeld experiment provided the case study data. Two sets of weather elements to represent different climatic conditions and four different soil types were speciﬁed, where results indicated little inﬂuence on these speciﬁcations for the identiﬁcation of the most sensitive parameters. Soil conditions and management were found to affect the ranking of parameter sensitivities more strongly than weather conditions for the selected outputs. Parameters related to drainage were strongly inﬂuential for simulations of soil water dynamics, yield and biomass of wheat, runoff, and leaching from soil during individual and consecutive growing years. Wheat yield and biomass simulations were sensitive to the ‘ammonium immobilised fraction’ parameter that related to soil mineralization and immobilisation. Simulations of CO 2 release from the soil and soil nutrient pool changes were most sensitive to external nutrient inputs and the process of denitriﬁcation, mineralization, and decomposition. This study provides important evidence of which SPACSYS parameters require the most care in their speciﬁcation. Moving forward, this evidence can help direct efﬁcient sampling and lab analyses for increased accuracy of such parameters. Results provide a useful reference for model users on which parameters are most inﬂuential for different simulation goals, which in turn provides better informed decision making for farmers and government policy alike. (fertiliser dissolution rate), LLF (loss of litter), and RLF (loss of residue) that are each related to external nutrient inputs also played a critical role in the simulations of surface runoff losses. Results showed that SFD was an inﬂuential parameter for the simulations of NHL (NH4 leach); LLF was inﬂuential for the simulations of NDR and CDR (dissolved N and C loss); and RLF was inﬂuential for the simulations of NRR and CRR (residue N and C loss). External nutrient inputs can increase soil nutrient pools, and in turn, nutrient losses when surface runoff occurs.


Introduction
Process-based models for agricultural systems provide a widely used and efficient tool for understanding the complex interactions between soil water, carbon (C), nitrogen (N), phosphorus (P), and plant growth [1,2], where both production under various environmental conditions [3] and nutrient cycling [4] can be simulated. To run a simulation successfully, models need to be informed with parameters that accurately quantify individual processes, together with local soil conditions, weather, and management practices. However, parameters that might vary with environmental conditions and have varying, possibly complex distributions themselves [5] are often difficult to precisely characterise, especially if they are costly to measure, requiring long-term monitoring [6]. In turn, this Agriculture 2021, 11, 624 2 of 30 parameter variability, uncertainty, or quality directly affects the reliability of the model simulations [1,7]. In general, agricultural simulation models require a large amount of input data and parameters, many of which are difficult to source or collect. At the same time, it is common for only a few parameters to strongly influence (and maximize) the variability of model outputs, while the majority of parameters only provide a weak influence in this respect [8]. Therefore, understanding the likely influence of each parameter on a model's outputs is crucial, as this allows the number of parameters to be safely reduced without a worrying loss of model accuracy, while at the same time, focus can be placed on data collection for parameters that most strongly influence model performance.
To identify how parameters impact model outputs, sensitivity diagnostics can be calculated through running multiple model simulations while varying the parameters across specific ranges [9]. Sensitivity analysis for model parameterization and outputs have been applied to various cropping systems, soil types and climate conditions using a range of models. For example, Specka et al. [10] applied a parameter sensitivity analysis for modelling above-ground biomass with regard to a future model calibration and an improved understanding of model response patterns; Jabloun et al. [11] assessed the sensitivity of outputs (crop yield and N leaching) from a crop simulation model, where 128 parameters were assessed. A sensitivity analysis quantifies those parameters which are most influential on model outputs, and in this respect can guide the efforts towards improving their accuracy as well as the model's output accuracy [12].
The SPACSYS model is implemented in a modular approach and provides a field-scale and weather-driven dynamic process-based simulation model of water, C and N cycling between plants, soils, and microbes [13]. The model has been widely used to assess the impact of climate change [14][15][16], tillage [17], fertilizer application [16,18], and different cultivars [19] on agricultural systems in terms of crop yields, C and N budgets, soil physical properties, and soil water redistribution. However, the model requires over 200 parameters for the simulation of various processes and plant growth and development. It would therefore be informative and beneficial to assess which parameters the SPACSYS outputs are the most and the least sensitive to, so that there can be guidelines for its use when all 200 parameters are not readily available.
A sensitivity analysis has also been used in building and understanding the structure of agro-ecosystem models [20,21]. Distinctions are made between a local and global sensitivity analysis (LSA and GSA, respectively), where the former is also known as a 'one at a time' (OAT) approach, while the latter considers multiple parameters at the same time. Thus, an LSA does not consider any interactions between parameters, while a GSA does. Both GSA [22] and LSA [23][24][25][26] have merit, where in ideal situations both should be applied and their results objectively assessed and compared [27,28]. In studies where uncertainty and interactions in model parameters are minimal, then an LSA may suffice [27] over the inherently more complex, GSA. Specifically, Link et al. [27] concludes that an LSA may suffice for situations where the ranking of model parameters is of importance, while GSA should be used for a precise attribution of variance in model outputs. Of course, it should be noted that results and insights from literature are always dependent on the characteristics of the process-based models employed.
Given this study is the first to identify parameter input sensitivity for the SPACSYS model, we chose to present results in two stages: first, an LSA as presented here, and second, a GSA (which is prep.) for presentation elsewhere. This two-stage reporting approach is followed for computational reasons, together with differences in associated interpretations, visualisations, and comparisons of the LSA and GSA outputs. Given SPACSYS requires over 200 parameters, assessing and reporting the sensitivity of them all would be problematic, for both an LSA and a GSA alike. Further, different subsets of parameters would be chosen according to their likely importance in simulating different processes under different background conditions. Because of computational and associated interpretation constraints, an LSA can investigate the sensitivity of outputs to a greater number of input parameters and investigate the distribution of each input parameter more  , than that feasible with a GSA. A GSA for, say, only four parameters entails the  added investigation of six 2-parameter interactions, four 3-parameter interactions, and one  4-parameter interaction, which themselves vary on how the parameter input distributions  are described (likely done simply, with say low, medium, and high values only). Given these thorny design and implementation issues for a GSA, for this study, we only consider the effects of a single parameter change on outputs via an LSA. This allows for a greater number of parameters to be investigated coupled with more detailed descriptions of the parameter distributions than would be viable within a GSA.
For this study, we carried out an LSA on the parameters that control soil water redistribution, and C and N cycling with the SPACSYS model, using a winter wheat field experiment conducted at Rothamsted Research in Harpenden, UK as the case study. To investigate potential influences of climatic and soil conditions on the analysis, simulations with the same configurations were run with different climatic and soil conditions. The study objectives were: (1) to identify those parameters that have the maximum influence on the simulation outputs with respect to wheat production and environment risks under UK soils and climatic conditions, and (2) to examine and understand the relationships between parameter sensitivity and model outputs.

Study Site
Data from a winter wheat field experiment was used for this study. The experiment was conducted for three growing seasons from 2011 to 2014 in the experimental plots of "Exhaustion Land", one of the classical long-term experiments at Rothamsted Research in Harpenden (51 • 49 N, 0 • 21 W and 128 m a.s.l.). The winter wheat cultivar was Xi-19. The soil is classified as a Chromic Luvisol (FAO classification) with a silty clay loam texture topsoil. Daily maximum and minimum air temperatures and precipitation during the growing years are shown in Figure 1. Compared to the mean climatic conditions over the growing season between 1981 and 2010 (mean temperature is 8.7 • C, precipitation 611mm), the first growing season and the third growing season were warmer and wetter (9.4 • C and 770 mm, 10.1 • C and 814 mm, respectively), but the second season was cooler and dryer (7.9 • C and 672 mm, respectively). The total sunshine hours of the three growing seasons (1314, 1268, and 1363 h, respectively) were all higher than the 1981 to 2010 mean (1243 h). Data for the study period were downloaded from the electronic Rothamsted Archive (e-RA, http://www.era.rothamsted.ac.uk/, accessed on 30 June 2021). greater number of input parameters and investigate the distribution of each input parameter more intensively, than that feasible with a GSA. A GSA for, say, only four parameters entails the added investigation of six 2-parameter interactions, four 3-parameter interactions, and one 4-parameter interaction, which themselves vary on how the parameter input distributions are described (likely done simply, with say low, medium, and high values only). Given these thorny design and implementation issues for a GSA, for this study, we only consider the effects of a single parameter change on outputs via an LSA. This allows for a greater number of parameters to be investigated coupled with more detailed descriptions of the parameter distributions than would be viable within a GSA. For this study, we carried out an LSA on the parameters that control soil water redistribution, and C and N cycling with the SPACSYS model, using a winter wheat field experiment conducted at Rothamsted Research in Harpenden, UK as the case study. To investigate potential influences of climatic and soil conditions on the analysis, simulations with the same configurations were run with different climatic and soil conditions. The study objectives were: (1) to identify those parameters that have the maximum influence on the simulation outputs with respect to wheat production and environment risks under UK soils and climatic conditions, and (2) to examine and understand the relationships between parameter sensitivity and model outputs.

Study Site
Data from a winter wheat field experiment was used for this study. The experiment was conducted for three growing seasons from 2011 to 2014 in the experimental plots of "Exhaustion Land", one of the classical long-term experiments at Rothamsted Research in Harpenden (51°49′ N, 0°21′ W and 128 m a.s.l.). The winter wheat cultivar was Xi-19. The soil is classified as a Chromic Luvisol (FAO classification) with a silty clay loam texture topsoil. Daily maximum and minimum air temperatures and precipitation during the growing years are shown in Figure 1. Compared to the mean climatic conditions over the growing season between 1981 and 2010 (mean temperature is 8.7 °C, precipitation 611mm), the first growing season and the third growing season were warmer and wetter (9.4 °C and 770 mm, 10.1 °C and 814 mm, respectively), but the second season was cooler and dryer (7.9 °C and 672 mm, respectively). The total sunshine hours of the three growing seasons (1314, 1268, and 1363 h, respectively) were all higher than the 1981 to 2010 mean (1243 h). Data for the study period were downloaded from the electronic Rothamsted Archive (e-RA, http://www.era.rothamsted.ac.uk/, accessed on 30 June 2021).

Parameters and Simulated Outputs
A total of 61 model required parameters, describing both water redistribution and C and N processes, were chosen to investigate the responses of the system. To reduce complexity in the analysis, parameters that relate to crop growth and development were set as the usual default values, and thus not part of the 61 parameters investigated. The ranges (maximum to minimum) of the chosen parameters were ± 50% of the default values used in previous studies [29][30][31]. Each parameter was set to one of 100 different values to ensure sampling points are uniformly distributed across the probability distribution. A total number of 6100 SPACSYS simulations were run, each with a 3-years length, meaning that only one of 61 parameters (Table A1) was investigated in turn, whilst all others were set at their default. Details on the linkage of each parameter to each process can be found in various applied SPACSYS model studies [13,32] and also in subsequent, developmental studies [31].
For the simulation outputs, a total of 27 variables were considered. Variables include those for grain yield and above-ground dry matter biomass of winter wheat, soil water dynamics, losses (runoff, leaching, and release) from soil and nutrient pool sizes (Table A2).

Sensitivity Analysis and Diagnostics
As indicated above, only the effect of a single parameter change on the simulated outputs was considered-i.e., an LSA was followed. We analysed sensitivity to the 61 parameters by ranking the values of the root mean square deviation (RMSD) of simulated outputs as follows: where D represents the simulated value using the default parameter values  and S i represents the simulated values as defined above, at the ith step of parameter changes and n is equal to 100 (the number of simulations). In order to ensure that all results from the SPACSYS runs, with each of 61 parameters varying in turn, can be assessed and objectively compered to each other, we used the normalized RMSD (NRMSD) as follows: where S max is the maximum value of all the simulated outputs for the whole set (100 simulations) of a single parameter over the entire simulated period. The NRMSD diagnostic combines RMSD values across multiple output variables to produce an overall measure of model sensitivity, while normalising for the different scales of the outputs. It can be interpreted as a fraction of the maximum values. The higher the NRMSD, the more sensitive the output variables are to the given parameter specification. A parameter that yields an NRMSD equal to 0% means it contributes very little on the simulated outputs and can be taken as a fixed parameter [33]. Simulations are considered sensitive to parameters yielding NRMSD values > 1%. In order to capture the effects of different weather conditions on parameter sensitivity, we calculated the NRMSD in each individual growing year (2011/12, 12/2013, and 13/2014) and consecutive years (2011/13, 12/2014, and 11/2014). Further clarity was provided with averaged NRMSDs for the individual growing years and consecutive years: where NRMSD IN and NRMSD CD represent the averaged NRMSD for the individual and the consecutive years, respectively. NRMSD 1 , NRMSD 2 , NRMSD 3 , NRMSD 12 , NRMSD 23 ,

Simulated Climate and Soil Data
In order to assess the impact of climate and soil type on the results of the sensitivity analysis, four additional simulation runs were built with the same configuration as the main simulation experiment described-replacing the weather over the simulation period with another set collected for South West England (50 • 46 10 N, 3 • 54 05 W) and replacing soil properties with (a) a sandy soil with 10% clay, 70% sand, and 20% silt content in the top 10 cm layer, with 12% clay, 72% sand, and 16% silt content in the 10-20 cm layer and with 15% clay, 75% sand, and 10% silt content in the 20-30 cm layer; (b) a loam soil with 20% clay, 40% sand, and 40% silt content in the top 10 cm layer, with 22% clay, 42% sand, and 36% silt content in the 10-20 cm layer and with 25% clay, 45% sand, and 30% silt content in the 20-30 cm layer; and (c) a clay soil with 80% clay, 10% sand, and 10% silt content in the top 10 cm layer, with 76% clay, 12% sand, and 12% silt content in the 10-20 cm layer and with 70% clay, 15% sand, and 15% silt content in the 20-30 cm layer. Additional simulations were configured as described in Section 2.2, where the use of default parameter values was considered reasonable. The influences of the parameters on the chosen output variable under alternative weather conditions are shown in Figure 3 and with the three different soil types in Figure 4 (sandy), Figure 5 (loam), and Figure 6 (clay). The results from these additional climate and soil simulations were hardly different to those from the main simulation experiment (Figure 2), where the sensitive parameters with NRMSD greater than 1% were almost identical. However, the number of fixed parameters reduced, especially in the sandy soil where none of the 61 parameters could be taken as fixed.

Results
As the most sensitive parameters had similar trends across all the simulations (main and additional), regardless of climate or soil type, we now focus our attention only on the relationships between the parameters and the output variables from the main simulation experiment. Here the greatest sensitivity to parameter specification was observed for simulated outputs grouped as: (1) soil water dynamics, (2) crop dry matter accumulation, and (3) N and C losses and pools. The influences of the parameters on the chosen output variable under alternative weather conditions are shown in Figure 3 and with the three different soil types in Figure  4 (sandy), Figure 5 (loam), and Figure 6 (clay). The results from these additional climate and soil simulations were hardly different to those from the main simulation experiment (Figure 2), where the sensitive parameters with NRMSD greater than 1% were almost identical. However, the number of fixed parameters reduced, especially in the sandy soil where none of the 61 parameters could be taken as fixed.      As the most sensitive parameters had similar trends across all the simulations (main and additional), regardless of climate or soil type, we now focus our attention only on the relationships between the parameters and the output variables from the main simulation experiment. Here the greatest sensitivity to parameter specification was observed for simulated outputs grouped as: (1) soil water dynamics, (2) crop dry matter accumulation, and (3) N and C losses and pools.

Soil Water Dynamics
The drainpipe level (DPL) parameter had the greatest influence on the simulations for ground water flow (GWF) and water content change (WSC) with NRMSD values > 7% (Figure 2), especially when the specified value of DPL was from the surface to the 2 m soil depth (Figure 7b,e,h). Here GWF increased and then plateaued at around 2 m, while WSC decreased moving from positive to negative values. For DPL values between 0 and −2 m, surface runoff (SRO) decreased from about 570 mm/year to values of zero. Likewise, other parameters relating to the drainage system, such as distance between drainpipes (DBD) and minimum roughness length (MRL) were similarly influential with NRMSD values > 1% and > 3% ( Figure 2) for GWF, WSC, and SRO outputs (Figure 7). Here GWF, WSC, and SRO outputs all decreased as the MRL parameter increased from its minimum (0.001 m) to maximum specification (1 m) ( Table A1 in Appendix A), while only GWF and SRO outputs decreased as the DBD parameter increased from its minimum (2 m) to maximum specification (100 m) ( Table A1). The influence of minimum hydraulic conductivity (MHC), runoff first order rate coefficient (RFO) and maximum surface storage (MSS) parameters on GWF, WSC, and SRO outputs are given in Figure A1 in Appendix A.
The drainpipe level (DPL) parameter had the greatest influence on the simulations for ground water flow (GWF) and water content change (WSC) with NRMSD values > 7% (Figure 2), especially when the specified value of DPL was from the surface to the 2 m soil depth (Figure 7b,e,h). Here GWF increased and then plateaued at around 2 m, while WSC decreased moving from positive to negative values. For DPL values between 0 and −2 m, surface runoff (SRO) decreased from about 570 mm/year to values of zero. Likewise, other parameters relating to the drainage system, such as distance between drainpipes (DBD) and minimum roughness length (MRL) were similarly influential with NRMSD values > 1% and > 3% ( Figure 2) for GWF, WSC, and SRO outputs (Figure 7). Here GWF, WSC, and SRO outputs all decreased as the MRL parameter increased from its minimum (0.001 m) to maximum specification (1 m) ( Table A1 in Appendix A), while only GWF and SRO outputs decreased as the DBD parameter increased from its minimum (2 m) to maximum specification (100 m) ( Table A1). The influence of minimum hydraulic conductivity (MHC), runoff first order rate coefficient (RFO) and maximum surface storage (MSS) parameters on GWF, WSC, and SRO outputs are given in Figure A1 in Appendix A.

Dry Matter of Winter Wheat
The simulation of dry matter of leaves, stems, and grains for winter wheat (LDM, SDM, and GDM) are particularly sensitive to parameters DBD, MRL, DPL, and AIF (ammonium immobilisation fraction) following their NRMSD values ( Figure 2). DBD was the most sensitive parameter (NRMSD > 3%) for LDM, SDM, and GDM in the individual years, with increases in dry matter simulation only for DBD < 40 m. However, SDM and GDM were most sensitive to AIF with NRMSD > 4%, where dry matter simulation increased slightly across the full range of the AIF values specified (i.e., 0 to 1). Increases in dry matter were found for DPL > −2 m, while a stepped but small increase in dry matter was observed for MRL across its minimum to maximum values specified ( Figure A2 in Appendix A).

Losses with Surface Runoff
The simulations for nitrate (NO3) loss with surface runoff (NOR) were also highly sensitive to DPL, with NRMSD values of 1.38% and 1.73% in the individual and pairs of consecutive years, respectively ( Figure 2). For DPL > −2 m, NOR displayed an increasing trend ( Figure 8a). As would be expected, simulations for NOR were also sensitive to the fertiliser dissolution rate (SFD), where NOR simulations displayed a non-linear increase for SFD < 0.2. (Figure 8b). With NRMSD values > 0%, the simulations of dissolved N loss (NDR) and C loss (CDR) with surface runoff were sensitive to the loss of litter parameter (LLF) (Figure 8c), and the simulations of residue N loss (NRR) and residue C loss (CRR) were sensitive to the loss of residue parameter (RLF) (Figure 8d).

Nitrogen and Carbon Leaching
With the specification of parameters DBD, MRL, DPL, MHC, and ESP (the empirical scale in pore shape), each had a clear influence on the simulations of leaching for NO 3 (NOL), NH 4 (NHL), N (NDL), and C (CDL) in the individual and pairs of consecutive years ( Figure 2). The same variable simulations were also sensitive to parameters AIF, LLF, RLF, DFF (a transferring fraction of decomposed fresh litter to dissolved organic matter), DFL (a transferring fraction of litter to dissolved organic matter), CWF (coefficient in the water function for decomposition), and DPR (potential decomposition rate of dissolved organic matter). DBD was a highly influential parameter for NOL and NHL with NRMSDs > 3% whilst DPL was highly influential for NDL and CDL with NRMSDs > 4%. Relationships between output simulations of NOL, NHL, NDL, and CDL and parameters DBD, MRL, DPL, AIF, MHC, ESP, DFF, DFL, CWF, DPR, LLF, and RLF are given in Figures A3 and A4 in Appendix A.

Gas Emissions
Simulations of the N 2 O emission rate (N2O) were most sensitive to changes in the maximum autotrophic nitrification rate (MAN), secondly, the NO production fraction from nitrification (PFF-NO), and thirdly, the Q10 temperature coefficient for denitrification (Q1D) (Figure 2). Specification of the latter two parameters also strongly influenced the NO emission rate (NOE) and N 2 emission rate (N2E) (Figure 2). In addition, N2O simulations were found to be relatively sensitive to DBD, MRL, and DPL ( Figure 2). Figure 9 displays the relationships between the nitrogenous gas emissions outputs (NOE, N2O, and N2E) and a range of values for parameters PFF-NO and Q1D. For the range of PFF-NO specified, the simulation of NOE, N2O, and N2E were erratic with no clear pattern. Further plots for NOE only are given in Figure A5 in Appendix A.
Parameters CWF, Q1M (Q10 temperature coefficient for mineralization), and ASF (assimilation factor), in general, provided the highest NRMSDs ( Figure 2) with respect to influencing CO 2 emissions from various soil organic C pools (DRE: dissolved release; FLR: fresh litter release; HRE: humus release; MRE: microbial release). However, the simulation of HRE was most sensitive to the humus potential decomposition rate (HPD) (NRMSD > 15%) (Figure 2). Figure 10 displays the relationships between CO 2 emissions DRE, FLR, HRE, MRE with respect to the top-ranking sensitive parameters CWF, Q1M, ASF over the simulation period. Almost all CO 2 emission outputs deceased as CWF, Q1M, and ASF parameters were increased through their minimum/maximums set (Table A1); the exception being MRE with ASF. Further CO 2 emission plots depicting their sensitivity to key parameters are given in Figures A6-A8 in Appendix A.

Changes of Soil C and N Pools
Although simulations of soil N and C pools in humus (NHP and CHP, respectively) were most sensitive to HPD (Figure 2), no discernible trend in this sensitivity was observed ( Figure A9a in Appendix A). Simulations of N and C dissolved (NDP and CDP, respectively) were most sensitive to CWF, DFL (dissolved fraction in litter), HFL (humus fraction from fresh litter), and DBD (Figure 2), where in this case, clear increasing trends in this sensitivity were observed, some plateauing ( Figure A9b,e for CWF and DBD, respectively), some not ( Figure A9c,d for DFL and HFL, respectively). Further soil N and C pools plots depicting their sensitivity to key parameters are given in Figures A10 and A11 in Appendix A.

Changes of Soil C and N Pools
Although simulations of soil N and C pools in humus (NHP and CHP, respectively) were most sensitive to HPD (Figure 2), no discernible trend in this sensitivity was observed ( Figure A9a in Appendix A). Simulations of N and C dissolved (NDP and CDP, respectively) were most sensitive to CWF, DFL (dissolved fraction in litter), HFL (humus fraction from fresh litter), and DBD (Figure 2), where in this case, clear increasing trends in this sensitivity were observed, some plateauing ( Figure A9b,e for CWF and DBD, respectively), some not ( Figure A9c,d for DFL and HFL, respectively). Further soil N and C pools plots depicting their sensitivity to key parameters are given in Figures A10 and A11 in Appendix A.

Discussion
Parameters describing drainage implementation (DBD and DPL) play a critical role in the performance of the agriculture system reflecting changes in the weather and soil condition. For this study, simulations for soil water dynamics, crop dry matter accumulation and N and C losses, were all found to be highly sensitive to the specification of DBD and DPL. Our results corroborate with Ballantine and Tanner [34], and Ritzema [35], where both studies demonstrated how the drainage system was an important externality

Discussion
Parameters describing drainage implementation (DBD and DPL) play a critical role in the performance of the agriculture system reflecting changes in the weather and soil condition. For this study, simulations for soil water dynamics, crop dry matter accumulation and N and C losses, were all found to be highly sensitive to the specification of DBD and DPL. Our results corroborate with Ballantine and Tanner [34], and Ritzema [35], where both studies demonstrated how the drainage system was an important externality for agricultural production to prevent waterlogging and salinization of the soil in arid and semi-arid regions and strategically control drainage of excess water from the soil profile. It has also been reported that water use efficiency for wheat with a controlled drainage with a varying depth was 40% higher than that with the conventional drainage, depending on the crop stage [36]. Tomic et al. [37] showed that a DBD value of 25 m provided the highest wheat yield among seven different drainage systems, which is in agreement with our result in that DBD is one of the most sensitive parameters for the dry matter of winter wheat. However, this might not be transferrable to other crops, e.g., drainage control was found to have no significant influence on soybean yield [38]. Designing an appropriate drainage system in humid areas can reduce N loss from an agricultural field [39]. By raising the water level in the soil profile, which is equivalent to increasing DPL in our study here, controlled drainage has the potential to increase denitrification in the anaerobic zone [40]. Design of subsurface drainpipes involves the determination of depth (DPL), spacing (DBD) and pipe diameter (DPD) [41]. As a parameter of a drainage system, however, DPD had no influence on any of the 27 selected output variables (Figure 2), which implied the effects of vertical water flow were stronger than horizontal water flow towards the pipes in our study's wheat field.
Our analysis suggests that simulations were sensitive to changes in MRL (minimum roughness length), where it is known to strongly affect water and heat fluxes on the soil surface [42,43]. Wang et al. [44] have shown that an accurate assessment of soil moisture might be problematic if MRL is inappropriately parameterised. However, our sensitivity analysis showed that SRI, HSG, and HCV (all related to evaporation and transpiration processes, see Table A1) had no influence on any of the 27 selected outputs.

Sensitive Parameters for Water Dynamics
Soil water dynamics is a key component for nutrient cycling and crop growth and development. Our study showed that soil water dynamics were sensitive to MHC (minimum hydraulic conductivity), DPL, DBD, and MRL (Figures 2 and A1). This is not surprising given that: (1) MHC is a key parameter to control the soil water infiltration process where water moves downward from the surface [45] and (2) soil hydraulic conductivity determines surface energy and water fluxes and then soil water content [46][47][48].

Sensitive Parameters for Yield and Biomass
For simplicity, we did not include in our sensitivity analysis, the plant biological parameters. Therefore, the identified sensitivities influence biomass accumulation and the grain yield of winter wheat indirectly via other processes. Here DBD, MRL, DPL, and AIF are all important parameters for yield and aboveground biomass simulation (where DBD, MRL, and DPL also influence water dynamics in the soil, as discussed above). Unarguably, soil water content itself plays an important role in yield and aboveground dry matter accumulation under non-irrigation conditions [49,50]. It regulates N cycling and its availability to crops is important, as this in turn, affects crop growth and yield formation [51][52][53]. AIF that can be immobilised when the immobilisation process occurs, impacts crop growth and yield through changing mineralisation, and in sequence, the status of soil N and C content [1,54].

Sensitive Parameters for Losses from Soil
The drainage parameters DPL and DBD also had significant influence on C and N losses through surface runoff and leaching from soils. In addition, SFD (fertiliser dissolution rate), LLF (loss of litter), and RLF (loss of residue) that are each related to external nutrient inputs also played a critical role in the simulations of surface runoff losses. Results showed that SFD was an influential parameter for the simulations of NHL (NH4 leach); LLF was influential for the simulations of NDR and CDR (dissolved N and C loss); and RLF was influential for the simulations of NRR and CRR (residue N and C loss). External nutrient inputs can increase soil nutrient pools, and in turn, nutrient losses when surface runoff occurs.
Clearly, parameters related to nitrification and denitrification processes have a great influence on nitrogenous gas emissions. A higher nitrification rate increases N availability for the denitrification process that is enhanced by Q1D (Q10 temperature coefficient for denitrification). Specification of MAN (maximum autotrophic nitrification rate) had the biggest influence on NOE (NO emission rate), and both PFF-NO (NO production fraction from nitrification) and Q1D influenced NOE, N2O (N 2 O emission rate), and N2E (N 2 emission rate). Conversely, simulations for soil CO2 release were sensitive to a different set of parameters that control nutrient cycling. Here, the specification for Q1M (Q10 temperature coefficient for mineralization), CWF (coefficient in the water function for decomposition), and ASF (assimilation factor) all had a significant influence on DRE (dissolved release), FLR (fresh litter release), HRE (humus release), and MRE (microbial release) (i.e., outputs variables that generate from decomposition from various soil C pools). Our results are supported by previous studies [55][56][57]. The soil organic matter decomposition rate depends on substrates [58,59], a potential decomposition rate and other abiotic and biotic factors including soil temperature (related to Q1M), soil moisture (i.e., CWF), and microbial activities (i.e., ASF) [60].

Conclusions
The local sensitivity analysis for the SPACSYS model using data from a winter wheat field experiment found that parameters related to drainage not only affected drainage fluxes but also grain yield and aboveground dry matter accumulation of winter wheat, and losses of water, C and N from soil under given weather and soil conditions. In addition, parameters related to mineralization and immobilisation processes influenced crop yield and biomass accumulation significantly. Further, parameters that control the nitrification and denitrification processes had a great influence on nitrogenous gas emissions. All results were found to be largely insensitive to different climatic and soil conditions, where identified parameter sensitivity remained relatively the same.
Overall, this study has provided evidence of which parameters require the most care in their specification and in turn, which parameters need to be determined accurately. This evidence can lead to efficient and cost-effective sampling and lab analyses for such parameters. Results provide a reference for model users on which parameters are key for different simulation goals. However, although this study has provided some important advances in model understanding and use, we have only considered the influence of single parameter change at a time on the model outputs. Interaction effects through a global sensitivity analysis were not assessed, and as such, are a focus of current and complementary research.    Figure A1. The relationships between ground water flow (GWF), water content change (WSC) and surface runoff (SRO) with changes in parameters with NRMSE > 1% (minimum hydraulic conductivity (MHC), ammonium immobilised fraction (AIF), specific fertiliser dissolution rate (SFD) and maximum autotrophic nitrification rate (MAN)).  Agriculture 2021, 11, x FOR PEER REVIEW 24 of 33 DBD: distance between drainpipes; MRL: minimum roughness length; DPL: drain pipe level; AIF: ammonium immobilised fraction; MHC: minimum hydraulic conductivity; and ESP: empirical scale in pore shape. Figure A4. The relationships between leached losses of N and C and the parameters with NRMSEs higher than 1% (CDL: C dissolve leach; NDL: N dissolve leach; DFF: dissolved fraction from fresh litter; DFL: dissolved fraction in litter; CWF: coefficient in water function in the process of organic matter decomposition; DPR: potential decomposition rate of dissolved organic matter; LLF: litter unit loss fraction; RLF: residue unit loss fraction). Figure A4. The relationships between leached losses of N and C and the parameters with NRMSEs higher than 1% (CDL: C dissolve leach; NDL: N dissolve leach; DFF: dissolved fraction from fresh litter; DFL: dissolved fraction in litter; CWF: coefficient in water function in the process of organic matter decomposition; DPR: potential decomposition rate of dissolved organic matter; LLF: litter unit loss fraction; RLF: residue unit loss fraction). Agriculture 2021, 11, x FOR PEER REVIEW 25 of 33 Figure A5. The relationships between NOE and sensitive parameters (DBD: distance between drainpipes; MRL: minimum roughness length; DPL: drain pipe level; AIF: ammonium immobilised fraction; MAN: maximum autotrophic nitrification rate; SCF: soil cover fraction prevent infiltration; CWF: coefficient in water function in the process of organic matter decomposition; BTD: base temp to unity denitrification; BTN: base temp to unity nitrification). Figure A5. The relationships between NOE and sensitive parameters (DBD: distance between drainpipes; MRL: minimum roughness length; DPL: drain pipe level; AIF: ammonium immobilised fraction; MAN: maximum autotrophic nitrification rate; SCF: soil cover fraction prevent infiltration; CWF: coefficient in water function in the process of organic matter decomposition; BTD: base temp to unity denitrification; BTN: base temp to unity nitrification).
Agriculture 2021, 11, x FOR PEER REVIEW 26 of 33 Figure A6. The relationships between carbon dissolved release (DRE) and sensitive parameters (DBD: distance between drainpipes; MRL: minimum roughness length; DPL: drain pipe level; MHC: minimum hydraulic conductivity; HPD: humus potential decomposition rate; LLF: litter unit loss fraction; OAW: optimal available water content; HFD: humus fraction from fresh litter; DFF: dissolved fraction from fresh litter; DFL: dissolved fraction in litter; RLT: residue to litter transfer fraction; FLD: fresh litter potential decomposition rate; RLF: residue unit loss fraction; DPR: dissolved potential decomposition rate; BTM: base temp to unity mineralization). Figure A6. The relationships between carbon dissolved release (DRE) and sensitive parameters (DBD: distance between drainpipes; MRL: minimum roughness length; DPL: drain pipe level; MHC: minimum hydraulic conductivity; HPD: humus potential decomposition rate; LLF: litter unit loss fraction; OAW: optimal available water content; HFD: humus fraction from fresh litter; DFF: dissolved fraction from fresh litter; DFL: dissolved fraction in litter; RLT: residue to litter transfer fraction; FLD: fresh litter potential decomposition rate; RLF: residue unit loss fraction; DPR: dissolved potential decomposition rate; BTM: base temp to unity mineralization).  Agriculture 2021, 11, x FOR PEER REVIEW 28 of 33 Figure A8. The relationships between microbial release (MRE) and sensitive parameters (DBD: distance between drainpipes; MRL: minimum roughness length; DPL: drain pipe level; MHC: minimum hydraulic conductivity; HPD: humus potential decomposition rate; HFL: partitioning fraction to humus from decomposed fresh litter; LLF: litter unit loss fraction; OAW: optimal available water content; HFD: partitioning fraction to humus from decomposed dissolved organic matter; DFF: dissolved fraction from fresh litter; DFL: dissolved fraction in litter; RLT: residue to litter transfer fraction; FLD: potential decomposition rate of fresh litter; RLF: Residue unit loss fraction; DPR: potential decomposition rate of dissolved organic matter; BTM: base temp to unity mineralization; MMR: microbial maintenance respiration rate). Figure A8. The relationships between microbial release (MRE) and sensitive parameters (DBD: distance between drainpipes; MRL: minimum roughness length; DPL: drain pipe level; MHC: minimum hydraulic conductivity; HPD: humus potential decomposition rate; HFL: partitioning fraction to humus from decomposed fresh litter; LLF: litter unit loss fraction; OAW: optimal available water content; HFD: partitioning fraction to humus from decomposed dissolved organic matter; DFF: dissolved fraction from fresh litter; DFL: dissolved fraction in litter; RLT: residue to litter transfer fraction; FLD: potential decomposition rate of fresh litter; RLF: Residue unit loss fraction; DPR: potential decomposition rate of dissolved organic matter; BTM: base temp to unity mineralization; MMR: microbial maintenance respiration rate). Agriculture 2021, 11, x FOR PEER REVIEW 29 of 33 Figure A9. Relationships between changes in soil organic C and N pools and the parameters they are most sensitive to over the simulation period. NHP: N humus pool; CHP: C humus pool; NDP: N dissolve pool; CDP: C dissolve pool; HPD: humus potential decomposition rate; CWF: coefficient in water function; DFL: dissolved fraction in litter; HFL: humus fraction from fresh litter; DBD: distance between drainpipes. Figure A10. The relationships between nitrogen dissolved pool (NDP) and sensitive parameters (MRL: minimum roughness length; DPL: drain pipe level; MHC: minimum hydraulic conductivity; HPD: humus potential decomposition rate; Q1M: Q10 value for mineralization ;ASF: assimilation factor; LLF: litter unit loss fraction; OAW: optimal available water content; DFF: dissolved fraction from fresh litter; FLD: potential decomposition rate of fresh litter; RLF: residue unit loss fraction; DPR: potential decomposition rate of dissolved organic matter; BTM: base temp to unity mineralization). Figure A9. Relationships between changes in soil organic C and N pools and the parameters they are most sensitive to over the simulation period. NHP: N humus pool; CHP: C humus pool; NDP: N dissolve pool; CDP: C dissolve pool; HPD: humus potential decomposition rate; CWF: coefficient in water function; DFL: dissolved fraction in litter; HFL: humus fraction from fresh litter; DBD: distance between drainpipes.
Agriculture 2021, 11, x FOR PEER REVIEW 29 of 33 Figure A9. Relationships between changes in soil organic C and N pools and the parameters they are most sensitive to over the simulation period. NHP: N humus pool; CHP: C humus pool; NDP: N dissolve pool; CDP: C dissolve pool; HPD: humus potential decomposition rate; CWF: coefficient in water function; DFL: dissolved fraction in litter; HFL: humus fraction from fresh litter; DBD: distance between drainpipes. Figure A10. The relationships between nitrogen dissolved pool (NDP) and sensitive parameters (MRL: minimum roughness length; DPL: drain pipe level; MHC: minimum hydraulic conductivity; HPD: humus potential decomposition rate; Q1M: Q10 value for mineralization ;ASF: assimilation factor; LLF: litter unit loss fraction; OAW: optimal available water content; DFF: dissolved fraction from fresh litter; FLD: potential decomposition rate of fresh litter; RLF: residue unit loss fraction; DPR: potential decomposition rate of dissolved organic matter; BTM: base temp to unity mineralization). Figure A10. The relationships between nitrogen dissolved pool (NDP) and sensitive parameters (MRL: minimum roughness length; DPL: drain pipe level; MHC: minimum hydraulic conductivity; HPD: humus potential decomposition rate; Q1M: Q10 value for mineralization; ASF: assimilation factor; LLF: litter unit loss fraction; OAW: optimal available water content; DFF: dissolved fraction from fresh litter; FLD: potential decomposition rate of fresh litter; RLF: residue unit loss fraction; DPR: potential decomposition rate of dissolved organic matter; BTM: base temp to unity mineralization).
Agriculture 2021, 11, x FOR PEER REVIEW 30 of 33 Figure A11. The relationships between carbon dissolved pool (CDP) simulations and sensitive parameters (MRL: minimum roughness length; DPL: drain pipe level; MHC: minimum hydraulic conductivity; CWF: coefficient in water function in the process of organic matter decomposition; Q1M: Q10 value for mineralization; ASF: assimilation factor; LLF: litter unit loss fraction; DFF: dissolved fraction from fresh litter; RLT: residue to litter transfer fraction; FLD: potential decomposition rate of fresh litter; RLF: residue unit loss fraction; DPR: potential decomposition rate of dissolved organic matter; BTM: base temp to unity mineralization). Figure A11. The relationships between carbon dissolved pool (CDP) simulations and sensitive parameters (MRL: minimum roughness length; DPL: drain pipe level; MHC: minimum hydraulic conductivity; CWF: coefficient in water function in the process of organic matter decomposition; Q1M: Q10 value for mineralization; ASF: assimilation factor; LLF: litter unit loss fraction; DFF: dissolved fraction from fresh litter; RLT: residue to litter transfer fraction; FLD: potential decomposition rate of fresh litter; RLF: residue unit loss fraction; DPR: potential decomposition rate of dissolved organic matter; BTM: base temp to unity mineralization).